ICO Thinning Example: UAS-SfM data from PSINF Ute Creek area
Wade Tinkham & George Woolsey
30 January, 2026
Example Data
The purpose of this analysis is to demonstrate how UAS-SfM aerial
point cloud data processed with cloud2trees
can be used to inform treatment planning and design. This analysis will
demonstrate how to process the spatial tree list following the methods
defined in Churchill et
al. (2016) for ICO thinning. The objective is to demonstrate how the
ICO method can be implemented to provide a cut-retention list based on
target stand average density targets for clump proportions.
UAS data collection missions were completed on the Pike-San Isabel National Forest (PSINF) in October 2025 as part of a demonstration to CFRI. Concurrent field data was collected by CFRI crews which will be used in later data validation efforts.
UAS Collection
- 130 ha flown in 1 day
- ~5,000 RGB images captured at 2.5 cm resolution
- Geo-rectified, wall-to-wall imagery
- Structure-from-Motion (SfM) point cloud with segmented trees
- Time commitment:
- 1.5 field days (7 people; but really 3 people since 4 were primarily observers)
- 1 day of processing (1 person)
- Using a ~$30K drone flown at 600 ft
- Processing time will increase with a smaller drone and a 400 ft flying altitude
Field Collection
- 26 one‑tenth–acre plots
- Approximately 1 plot per 10–15 acres, depending on stand
density
- Standard forestry metrics with stem‑mapped overstory trees and
saplings
- Time commitment:
- ~7 field days
- 3–4 people
load the standard libraries we use to do work
# bread-and-butter
library(tidyverse) # the tidyverse
library(viridis) # viridis colors
library(harrypotter) # hp colors
library(RColorBrewer) # brewer colors
library(scales) # work with number and plot scales
library(latex2exp)
# visualization
library(mapview) # interactive html maps
library(kableExtra) # tables
library(patchwork) # combine plots
library(ggnewscale) # new scale
library(ggrepel) # repel labels
# spatial analysis
library(terra) # raster
library(sf) # simple features
library(dbscan) # spatial clusteringnow, define a function to convert metric to imperial units
calc_imperial_units_fn <- function(df) {
df %>%
# convert to imperial units
dplyr::mutate(
dplyr::across(
.cols = tidyselect::ends_with("_cm")
, ~ .x * 0.394
, .names = "{.col}_in"
)
, dplyr::across(
.cols = tidyselect::ends_with("_m")
, ~ .x * 3.28
, .names = "{.col}_ft"
)
, dplyr::across(
.cols = tidyselect::ends_with("_m2_per_ha")
, ~ .x * 4.359
, .names = "{.col}_ftac"
)
, dplyr::across(
.cols = tidyselect::ends_with("_per_ha") & !tidyselect::ends_with("_m2_per_ha")
, ~ .x * 0.405
, .names = "{.col}_ac"
)
, dplyr::across(
.cols = tidyselect::ends_with("_area_ha")
, ~ .x * 2.471
, .names = "{.col}_ac"
)
, dplyr::across(
.cols = tidyselect::ends_with("_m2")
, ~ .x * 10.764
, .names = "{.col}_ft2"
)
) %>%
dplyr::rename_with(
.fn = function(x){dplyr::case_when(
stringr::str_ends(x,"_cm_in") ~ stringr::str_replace(x,"_cm_in","_in")
, stringr::str_ends(x,"_m_ft") ~ stringr::str_replace(x,"_m_ft","_ft")
, stringr::str_ends(x,"_m2_per_ha_ftac") ~ stringr::str_replace(x,"_m2_per_ha_ftac","_ft2_per_ac")
, stringr::str_ends(x,"_per_ha_ac") ~ stringr::str_replace(x,"_per_ha_ac","_per_ac")
, stringr::str_ends(x,"_area_ha_ac") ~ stringr::str_replace(x,"_area_ha_ac","_area_ac")
, stringr::str_ends(x,"_m2_ft2") ~ stringr::str_replace(x,"_m2_ft2","_ft2")
, TRUE ~ x
)}
)
}Load Data
load data
# where is the processed data from cloud2trees?
input_dir <- "../data/"
# tree top points
treetops_sf_with_dbh <- sf::st_read(file.path(input_dir, "West_Stand_final_detected_tree_points.gpkg"))## Reading layer `West_Stand_final_detected_tree_points' from data source
## `C:\Data\usfs\psinf_utecreek_ico_ex\data\West_Stand_final_detected_tree_points.gpkg'
## using driver `GPKG'
## Simple feature collection with 13964 features and 61 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 465870.1 ymin: 4355809 xmax: 466830.4 ymax: 4356255
## Projected CRS: WGS 84 / UTM zone 13N# what?
treetops_sf_with_dbh %>% dplyr::glimpse()## Rows: 13,964
## Columns: 62
## $ treeID <chr> "57060_466401.1_4355851.6", "27967_466607.1_…
## $ tree_height_m <dbl> 24.723, 24.350, 23.548, 22.951, 23.115, 22.8…
## $ crown_area_m2 <dbl> 41.8125, 32.8750, 20.5000, 18.5625, 53.0625,…
## $ fia_est_dbh_cm <dbl> 34.92759, 34.79764, 34.49266, 34.27861, 34.2…
## $ fia_est_dbh_cm_lower <dbl> 19.41092, 19.13628, 19.17295, 18.98428, 19.0…
## $ fia_est_dbh_cm_upper <dbl> 56.16955, 56.11775, 56.01807, 55.60062, 55.2…
## $ dbh_cm <dbl> 34.92759, 34.79764, 34.49266, 34.27861, 34.2…
## $ is_training_data <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
## $ dbh_m <dbl> 0.3492759, 0.3479764, 0.3449266, 0.3427861, …
## $ radius_m <dbl> 0.1746379, 0.1739882, 0.1724633, 0.1713930, …
## $ basal_area_m2 <dbl> 0.09581358, 0.09510194, 0.09344225, 0.092286…
## $ basal_area_ft2 <dbl> 1.0313374, 1.0236773, 1.0058124, 0.9933675, …
## $ ptcld_extracted_dbh_cm <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ ptcld_predicted_dbh_cm <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ tree_cbh_m <dbl> 8.388833, 9.573767, 8.006800, 8.593433, 7.82…
## $ is_training_cbh <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
## $ forest_type_group_code <chr> "220", "220", "220", "220", "220", "220", "2…
## $ forest_type_group <chr> "Ponderosa pine group", "Ponderosa pine grou…
## $ hardwood_softwood <chr> "Softwood", "Softwood", "Softwood", "Softwoo…
## $ comp_trees_per_ha <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_relative_tree_height <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_dist_to_nearest_m <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ max_crown_diam_height_m <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ is_training_hmd <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ landfire_stand_id <dbl> 2144, 1468, 1460, 1583, 1974, 1526, 2145, 21…
## $ crown_dia_m <dbl> 7.296391, 6.469757, 5.108954, 4.861534, 8.21…
## $ crown_length_m <dbl> 16.33417, 14.77623, 15.54120, 14.35757, 15.2…
## $ crown_volume_m3 <dbl> 455.3149, 323.8458, 212.3964, 177.6749, 540.…
## $ landfire_tree_kg_per_m3 <dbl> 0.1890150, 0.2128839, 0.2036069, 0.2991796, …
## $ landfire_stand_kg_per_m3 <dbl> 0.08, 0.11, 0.11, 0.17, 0.17, 0.17, 0.08, 0.…
## $ landfire_crown_biomass_kg <dbl> 86.06134, 68.94156, 43.24538, 53.15671, 154.…
## $ FIAnum <int> 122, 122, 122, 122, 122, 122, 122, 122, 122,…
## $ REGION <chr> "02", "02", "02", "02", "02", "02", "02", "0…
## $ FORESTNUMB <chr> "12", "12", "12", "12", "12", "12", "12", "1…
## $ DISTRICTNU <chr> "11", "11", "11", "11", "11", "11", "11", "1…
## $ spcd <int> 122, 122, 122, 122, 122, 122, 122, 122, 122,…
## $ voleq <chr> "200FW2W122", "200FW2W122", "200FW2W122", "2…
## $ TCFV <dbl> 31.2, 30.5, 28.9, 27.8, 28.0, 27.5, 27.4, 27…
## $ CFV_GRS <dbl> 28.9, 28.9, 26.0, 25.6, 26.0, 25.6, 25.6, 25…
## $ VStump <dbl> 0.9291809, 0.9232748, 0.9086702, 0.8982098, …
## $ VTip <dbl> 0.6295683, 0.5614567, 0.6773489, 0.6989350, …
## $ BFV_GRS <dbl> 150, 150, 130, 130, 130, 130, 130, 130, 120,…
## $ BFV_GRS_INTL <dbl> 165, 165, 140, 140, 140, 140, 140, 140, 130,…
## $ Poly_Name <chr> "West", "West", "West", "West", "West", "Wes…
## $ Desc_ <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ Poly <chr> "Boundary Adjusted", "Boundary Adjusted", "B…
## $ CN <chr> "c6916fad-74f1-4326-a788-f3c2eaf388b9", "c69…
## $ Area_MtSq <dbl> 69250.85, 69250.85, 69250.85, 69250.85, 6925…
## $ Area_Ac <dbl> 71.63097, 71.63097, 71.63097, 71.63097, 71.6…
## $ Area_Ha <dbl> 6.925085, 6.925085, 6.925085, 6.925085, 6.92…
## $ Perim_M <dbl> 1369.383, 1369.383, 1369.383, 1369.383, 1369…
## $ PerimL_M <dbl> 817.6369, 817.6369, 817.6369, 817.6369, 817.…
## $ Desc1 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ Shape_Leng <dbl> 2589.825, 2589.825, 2589.825, 2589.825, 2589…
## $ Shape_Area <dbl> 289880.3, 289880.3, 289880.3, 289880.3, 2898…
## $ size_class <chr> "Saw Timber (9-14.9 inch DBH)", "Saw Timber …
## $ dbscan_cluster <chr> "0", "6", "6", "6", "6", "6", "6", "6", "6",…
## $ clump_id <chr> "339", "6", "6", "6", "6", "6", "6", "6", "6…
## $ clump_n_trees <int> 1, 12235, 12235, 12235, 12235, 12235, 12235,…
## $ clump_class <chr> "Individual", ">15 trees", ">15 trees", ">15…
## $ largest_in_clump <chr> "Largest", "Largest", "Second Largest", "", …
## $ geom <POINT [m]> POINT (466401.1 4355852), POINT (46660…generate a stand boundary since we don’t have one at the moment and look at where we are in the US
stand_sf <-
treetops_sf_with_dbh %>%
sf::st_union() %>%
sf::st_concave_hull(ratio = 0.08) %>%
sf::st_buffer(2) %>%
sf::st_as_sf() %>%
dplyr::mutate(dummy=1) %>%
sf::st_as_sfc()
# where?
mapview::mapview(
stand_sf
, color = "black"
, lwd = 2
, alpha.regions = 0
, label = FALSE
, legend = FALSE
, popup = FALSE
, layer.name = "stand boundary"
)let’s make a base ggplot2 of the stand
plt_stand <- ggplot2::ggplot() +
ggplot2::geom_sf(data = stand_sf, color = "black", fill = NA) +
ggplot2::scale_x_continuous(expand = c(0, 0)) +
ggplot2::scale_y_continuous(expand = c(0, 0)) +
ggplot2::labs(
x = ""
, y = ""
) +
ggplot2::theme_void() +
ggplot2::theme(
legend.position = "top" # c(0.5,1)
, legend.direction = "horizontal"
, legend.margin = ggplot2::margin(0,0,0,0)
, legend.text = ggplot2::element_text(size = 8)
, legend.title = ggplot2::element_text(size = 8)
, legend.key = ggplot2::element_rect(fill = NA, color = NA)
# , plot.title = ggtext::element_markdown(size = 10, hjust = 0.5)
, plot.title = ggplot2::element_text(size = 10, hjust = 0.5, face = "bold")
, plot.subtitle = ggplot2::element_text(size = 8, hjust = 0.5, face = "italic")
)
# plt_standDefine Overstory and Clump Spacing
let’s define our DBH limit for classifying trees as “overstory” and we’ll define our distance in meters to separate tree clump groupings
See Churchill et al. (2016) Figure 4 (p.10) and Matonis and Binkley (2018) who “calculated coverage of mosaic-meadows (percentage of stand > 6 m from overstory trees)” (p. 124)
my_clump_dist_m <- 6
my_ostory_dbh_cm <- 5*2.54 # = inches*2.54plot the stand with some trees
plt_stand +
ggplot2::geom_sf(
data = treetops_sf_with_dbh %>%
sf::st_intersection(stand_sf) %>%
# overstory
dplyr::mutate(overstory = dbh_cm>=my_ostory_dbh_cm) %>%
dplyr::arrange(overstory)
, mapping = ggplot2::aes(color = overstory)
, size = 0.8, alpha = 0.7
) +
ggplot2::scale_color_manual(values=c("gray","navy"))
ggplot2::ggsave("../data/01_tree_tops.jpg", width = 7, height = 5)Define Clumps
create a general function to apply dbscan::dbscan to
point data of class sf
With respect to clump size groupings, Churchill et al. (2016) note that:
Proportions for clump sizes should be lumped into four or five bins for operational simplicity. We use 4 or 5 bins (Fig 5): individual trees, small clumps (2-4 trees), medium clumps (5-9 trees), and large clumps (10-20+ trees). Note that when instructed to leave a large clump (e.g. 10-20 trees), marking crews often have difficulty leaving the upper end of the size range (e.g. an 18, 19, or 20 tree clump). Thus adding a fifth bin for “super clumps” may be necessary (e.g. 15-20 trees or 20-25+ trees), especially if the upper size range of clumps is desired. (p. 12-13)
In their prescription worksheet Excel file Churchill et al. (2016) have the following group sizes: 1, 2-4 (Mean clump size: 3), 5-9 (7), 10-15 (12), 16-25 (20). To be able to cut clumps down to achieve the desired largest “16-25” group, we’ll add a “super” group with >25 trees. This same principle can be applied if you desire a clump size group larger than this (i.e. 25-35), whereby a group above this should be added so that the resulting prescription cuts the “super” group to ensure the proper sizing of the desired largest group.
# function to clump data that are sf points
st_clump_points <- function(
x # point data of class `sf`
, clump_dist_m = 6 # size (radius) of the epsilon neighborhood = maximum distance between points to add to clump
, clump_breaks = c(0,1,4,9,15,25,Inf) # where to break the clump size groups
, clump_labels = c("Individual","2-4 trees","5-9 trees","10-15 trees","16-25 trees",">25 trees") # what to name the clump size groups
) {
if(!inherits(x,"sf")){stop("`x` requires sf data")}
if(
is.na(clump_dist_m) ||
is.null(clump_dist_m) ||
!inherits(as.numeric(clump_dist_m),"numeric")
){
stop("`clump_dist_m` must be numeric and not missing")
}else{clump_dist_m <- as.numeric(clump_dist_m)[1]}
if(
any(is.na(clump_breaks)) ||
any(is.null(clump_breaks)) ||
!inherits(as.numeric(clump_breaks),"numeric")
){
stop("`clump_breaks` must be numeric and not missing")
}else{clump_breaks <- as.numeric(clump_breaks)}
# get points as x,y
point_clusters <- x %>%
dplyr::mutate(
X = sf::st_coordinates(.)[,1] %>% as.numeric()
, Y = sf::st_coordinates(.)[,2] %>% as.numeric()
) %>%
# get rid of columns we'll create
dplyr::select( -dplyr::any_of(c(
"hey_xxxxxxxxxx"
, "dbscan_cluster"
, "clump_id"
, "clump_n_trees"
, "clump_n_trees_grp"
, "clump_dist_m"
, "clump_n_trees_lb"
, "clump_n_trees_ub"
)))
#############################################################################
##### Identify clusters in each stem map plot #####
#############################################################################
### Place trees into clusters using an inter-tree distance of 6 m
my_dbscan_temp <- point_clusters %>%
sf::st_drop_geometry() %>%
dplyr::select(X,Y) %>%
dbscan::dbscan(eps = clump_dist_m, minPts = 2)
# my_dbscan_temp %>% str()
### append cluster ID to tree points
point_clusters$dbscan_cluster <- my_dbscan_temp$cluster
# point_clusters$cluster %>% summary()
# point_clusters %>% sf::st_drop_geometry() %>% dplyr::count(cluster) %>% dplyr::arrange(desc(n)) %>% dplyr::slice_head(n=11)
# clump_labels
if( !inherits(clump_labels,"character") ){
clump_labels <- NULL
}
### cluster metrics
point_clusters <- point_clusters %>%
dplyr::group_by(dbscan_cluster) %>%
dplyr::mutate(
# unique dbscan_cluster for individuals
clump_id = dplyr::case_when(
# cluster = 0 from dbscan::dbscan() are outliers (e.g. individual trees)
dbscan_cluster == 0 ~ max(my_dbscan_temp$cluster)+dplyr::row_number()
, T ~ dbscan_cluster
) %>%
factor()
) %>%
dplyr::group_by(clump_id) %>%
dplyr::mutate(
dbscan_cluster = factor(dbscan_cluster)
, clump_n_trees = dplyr::n()
, clump_n_trees_grp = base::cut(
clump_n_trees
, breaks = clump_breaks
, labels = NULL
)
) %>%
dplyr::ungroup() %>%
dplyr::mutate(clump_dist_m = clump_dist_m) %>%
dplyr::select(-c(dbscan_cluster))
### clump factor lb and ub
point_clusters <- point_clusters %>%
dplyr::bind_cols(
# get the cut upper and lower bound
tidyr::separate_wider_regex(
point_clusters
, cols = clump_n_trees_grp
, patterns = c(
"[\\(\\[]"
, clump_n_trees_lb = ".*"
, ","
, clump_n_trees_ub = ".*"
, "[\\)\\]]"
)
, cols_remove = FALSE
) %>%
# "Inf" string becomes numeric Inf here
dplyr::mutate(
# dplyr::across(
# c(clump_n_trees_lb, clump_n_trees_ub)
# , as.numeric
# )
# , clump_n_trees_lb = clump_n_trees_lb+1
# store as factor so we can use cut levels later
clump_n_trees_lb = factor(
as.numeric(clump_n_trees_lb)+1
, levels = as.character( clump_breaks[-length(clump_breaks)]+1 ) ## all but the last
# , ordered = T
)
, clump_n_trees_ub = factor(
clump_n_trees_ub
, levels = as.character(clump_breaks[-1]) ## all but the first
# , ordered = T
)
) %>%
dplyr::select(clump_n_trees_lb, clump_n_trees_ub)
)
### clump factor
if(!any(is.null(clump_labels))){
point_clusters <- point_clusters %>%
dplyr::mutate(
clump_n_trees_grp = clump_n_trees_grp %>%
forcats::fct_relabel(
~ clump_labels[match(.x, levels(clump_n_trees_grp))]
)
)
}
# ensure final clump factor is ordered ... nah, this messes up when levels aren't in data
# point_clusters <-
# point_clusters %>%
# dplyr::mutate(
# clump_n_trees_grp = clump_n_trees_grp %>%
# # order by the lower bound just to be sure
# forcats::fct_reorder(as.numeric(as.character(clump_n_trees_lb))) %>%
# # ordered so can make >,< comparisions
# ordered()
# )
# return
return(point_clusters)
}
# st_clump_points(
# treetops_sf_with_dbh %>% dplyr::slice_sample(n=1111)
# # treetops_clustered
# , clump_dist_m = 6
# , clump_breaks = c(0,1,9,Inf)
# , clump_labels = NULL
# ) %>%
# # dplyr::glimpse()
# sf::st_drop_geometry() %>%
# dplyr::count(clump_n_trees_grp,clump_n_trees_lb,clump_n_trees_ub)function to filter the tree points attached to stand bounds and apply
the st_clump_points function
# create function to pass a unit id and return list of trees with clump groupings
get_tree_clumps <- function(
trees_sf
, stand_sf
, ostory_ht_m = NA
, ostory_dbh_cm = NA
, clump_dist_m = NA
, clump_breaks = NA
, clump_labels = NA
){
if(
is.na(clump_dist_m) ||
is.null(clump_dist_m) ||
!inherits(as.numeric(clump_dist_m),"numeric")
){
stop("`clump_dist_m` must be numeric and not missing")
}else{clump_dist_m <- as.numeric(clump_dist_m)[1]}
if(
any(is.na(clump_breaks)) ||
any(is.null(clump_breaks)) ||
!inherits(as.numeric(clump_breaks),"numeric")
){
stop("`clump_breaks` must be numeric and not missing")
}else{clump_breaks <- as.numeric(clump_breaks)}
# check data
if(!inherits(stand_sf,"sf") && !inherits(stand_sf,"sfc")){stop("`stand_sf` requires polygon sf data")}
if(
!all( sf::st_is(stand_sf, c("POLYGON","MULTIPOLYGON")) )
){
stop("`stand_sf` requires POLYGON sf data")
}else{
stand_sf <- sf::st_union(stand_sf)
}
if(!inherits(trees_sf,"sf")){stop("`trees_sf` requires polygon sf data")}
if(
!all( sf::st_is(trees_sf, c("POINT","MULTIPOINT")) )
){
stop("`trees_sf` requires POINT sf data")
}
# check ostory definition
if(is.na(as.numeric(ostory_dbh_cm)) && is.na(as.numeric(ostory_ht_m))){
warning("`ostory_dbh_cm` and `ostory_ht_m` are not set...using `ostory_dbh_cm` = 12.7")
if(
!any(stringr::str_equal(names(trees_sf), "dbh_cm"))
){
stop("`trees_sf` requires 'dbh_cm' column")
}
ostory_dbh_cm <- 5*2.54
# filter data
ttops_temp <- trees_sf %>%
sf::st_intersection(stand_sf %>% sf::st_transform(sf::st_crs(trees_sf))) %>%
dplyr::filter(
dbh_cm>=as.numeric(ostory_dbh_cm)
)
}else if(!is.na(as.numeric(ostory_dbh_cm))){
if(
!any(stringr::str_equal(names(trees_sf), "dbh_cm"))
){
stop("`trees_sf` requires 'dbh_cm' column")
}
# filter data
ttops_temp <- trees_sf %>%
sf::st_intersection(stand_sf %>% sf::st_transform(sf::st_crs(trees_sf))) %>%
dplyr::filter(
dbh_cm>=as.numeric(ostory_dbh_cm)
)
}else{
if(
!any(stringr::str_equal(names(trees_sf), "tree_height_m"))
){
stop("`trees_sf` requires 'tree_height_m' column")
}
# filter data
ttops_temp <- trees_sf %>%
sf::st_intersection(stand_sf %>% sf::st_transform(sf::st_crs(trees_sf))) %>%
dplyr::filter(
tree_height_m>=as.numeric(ostory_ht_m)
)
}
### Place trees into clusters using an inter-tree distance
ttops_temp <- st_clump_points(
x = ttops_temp
, clump_dist_m = clump_dist_m
, clump_breaks = clump_breaks
, clump_labels = clump_labels
)
# get rid of columns we'll create
ttops_temp <-
ttops_temp %>%
dplyr::select( -dplyr::any_of(c(
"hey_xxxxxxxxxx"
, "clump_nn_dist_m"
, "clump_ostory_ht_m"
, "clump_ostory_dbh_cm"
)))
# add distance to nearest tree within clump
#### !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!GW: update to use kdtree:
#### !!!!!!!!!!!!!!!!!!! the `RANN` [package](https://jefferislab.github.io/RANN/)
#### !!!!!!!!!!!!!!!!!!! enables KD-tree processing using X and Y coordinates
#### !!!!!!!!!!!!!!!!!!! to build the tree which is a super-fast way to find the `x` number of
#### !!!!!!!!!!!!!!!!!!! near neighbors for each point in an input dataset (see `RANN::nn2()`).
ttops_temp <-
ttops_temp %>%
dplyr::group_by(clump_id) %>%
tidyr::nest() %>%
dplyr::mutate(
clump_nn_dist_m = purrr::map(data, function(x){
# get index of nearest neighbor
i = sf::st_nearest_feature(x)
# get dist
d = sf::st_distance(x, x[i,], by_element=TRUE) %>% as.numeric()
return(d)
})
) %>%
tidyr::unnest(cols = c(data, clump_nn_dist_m)) %>%
sf::st_set_geometry("geom") %>% # set it cuz it got lost
dplyr::ungroup() %>%
dplyr::mutate(
clump_dist_m = clump_dist_m
# , clump_id_duplicate = clump_id # can use this even after nesting data by clump_id
, clump_ostory_ht_m = ifelse(is.na(ostory_ht_m), as.numeric(NA), as.numeric(ostory_ht_m))
, clump_ostory_dbh_cm = ifelse(is.na(ostory_dbh_cm), as.numeric(NA), as.numeric(ostory_dbh_cm))
) %>%
dplyr::relocate(
tidyselect::starts_with("clump_")
, .after = dplyr::last_col()
)
# return
return(ttops_temp)
}call that get_tree_clumps() function which uses our
st_clump_points() function for overstory trees only
# call it
treetops_clustered <-
get_tree_clumps(
trees_sf = treetops_sf_with_dbh # %>% dplyr::slice_sample(n=2222)
, stand_sf = stand_sf
, ostory_dbh_cm = my_ostory_dbh_cm
, clump_dist_m = my_clump_dist_m
, clump_breaks = c(0,1,4,9,15,25,Inf) # where to break the clump size groups
, clump_labels = c("Individual","2-4 trees","5-9 trees","10-15 trees","16-25 trees",">25 trees") # what to name the clump size groups
)what?
treetops_clustered %>% dplyr::glimpse()## Rows: 8,808
## Columns: 68
## $ treeID <chr> "57060_466401.1_4355851.6", "27967_466607.1_…
## $ tree_height_m <dbl> 24.723, 24.350, 21.165, 21.191, 20.848, 20.9…
## $ crown_area_m2 <dbl> 41.8125, 32.8750, 19.0625, 31.6875, 29.4375,…
## $ fia_est_dbh_cm <dbl> 34.92759, 34.79764, 33.62946, 33.37739, 33.3…
## $ fia_est_dbh_cm_lower <dbl> 19.41092, 19.13628, 18.30126, 18.45933, 18.4…
## $ fia_est_dbh_cm_upper <dbl> 56.16955, 56.11775, 54.53841, 54.07352, 54.0…
## $ dbh_cm <dbl> 34.92759, 34.79764, 33.62946, 33.37739, 33.3…
## $ is_training_data <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
## $ dbh_m <dbl> 0.3492759, 0.3479764, 0.3362946, 0.3337739, …
## $ radius_m <dbl> 0.1746379, 0.1739882, 0.1681473, 0.1668870, …
## $ basal_area_m2 <dbl> 0.09581358, 0.09510194, 0.08882384, 0.087497…
## $ basal_area_ft2 <dbl> 1.0313374, 1.0236773, 0.9560998, 0.9418210, …
## $ ptcld_extracted_dbh_cm <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ ptcld_predicted_dbh_cm <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ tree_cbh_m <dbl> 8.388833, 9.573767, 8.747933, 9.090733, 8.29…
## $ is_training_cbh <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
## $ forest_type_group_code <chr> "220", "220", "220", "220", "220", "220", "2…
## $ forest_type_group <chr> "Ponderosa pine group", "Ponderosa pine grou…
## $ hardwood_softwood <chr> "Softwood", "Softwood", "Softwood", "Softwoo…
## $ comp_trees_per_ha <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_relative_tree_height <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_dist_to_nearest_m <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ max_crown_diam_height_m <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ is_training_hmd <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ landfire_stand_id <dbl> 2144, 1468, 1411, 1468, 1527, 1468, 1412, 14…
## $ crown_dia_m <dbl> 7.296391, 6.469757, 4.926574, 6.351833, 6.12…
## $ crown_length_m <dbl> 16.334166, 14.776234, 12.417068, 12.100267, …
## $ crown_volume_m3 <dbl> 455.31488, 323.84579, 157.80023, 255.61813, …
## $ landfire_tree_kg_per_m3 <dbl> 0.1890150, 0.2128839, 0.1841485, 0.2128839, …
## $ landfire_stand_kg_per_m3 <dbl> 0.08, 0.11, 0.11, 0.11, 0.12, 0.11, 0.12, 0.…
## $ landfire_crown_biomass_kg <dbl> 86.06134, 68.94156, 29.05867, 54.41699, 40.1…
## $ FIAnum <int> 122, 122, 122, 122, 122, 122, 122, 122, 122,…
## $ REGION <chr> "02", "02", "02", "02", "02", "02", "02", "0…
## $ FORESTNUMB <chr> "12", "12", "12", "12", "12", "12", "12", "1…
## $ DISTRICTNU <chr> "11", "11", "11", "11", "11", "11", "11", "1…
## $ spcd <int> 122, 122, 122, 122, 122, 122, 122, 122, 122,…
## $ voleq <chr> "200FW2W122", "200FW2W122", "200FW2W122", "2…
## $ TCFV <dbl> 31.2, 30.5, 25.0, 24.6, 24.2, 24.2, 24.1, 22…
## $ CFV_GRS <dbl> 28.9, 28.9, 21.5, 21.5, 21.5, 21.5, 21.5, 20…
## $ VStump <dbl> 0.9291809, 0.9232748, 0.8594608, 0.8452151, …
## $ VTip <dbl> 0.6295683, 0.5614567, 0.5501549, 0.5516371, …
## $ BFV_GRS <dbl> 150, 150, 110, 110, 110, 110, 110, 110, 110,…
## $ BFV_GRS_INTL <dbl> 165, 165, 125, 125, 125, 125, 125, 120, 120,…
## $ Poly_Name <chr> "West", "West", "West", "West", "West", "Wes…
## $ Desc_ <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ Poly <chr> "Boundary Adjusted", "Boundary Adjusted", "B…
## $ CN <chr> "c6916fad-74f1-4326-a788-f3c2eaf388b9", "c69…
## $ Area_MtSq <dbl> 69250.85, 69250.85, 69250.85, 69250.85, 6925…
## $ Area_Ac <dbl> 71.63097, 71.63097, 71.63097, 71.63097, 71.6…
## $ Area_Ha <dbl> 6.925085, 6.925085, 6.925085, 6.925085, 6.92…
## $ Perim_M <dbl> 1369.383, 1369.383, 1369.383, 1369.383, 1369…
## $ PerimL_M <dbl> 817.6369, 817.6369, 817.6369, 817.6369, 817.…
## $ Desc1 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ Shape_Leng <dbl> 2589.825, 2589.825, 2589.825, 2589.825, 2589…
## $ Shape_Area <dbl> 289880.3, 289880.3, 289880.3, 289880.3, 2898…
## $ size_class <chr> "Saw Timber (9-14.9 inch DBH)", "Saw Timber …
## $ X <dbl> 466401.1, 466607.1, 466600.6, 466608.1, 4666…
## $ Y <dbl> 4355852, 4356234, 4356244, 4356218, 4356208,…
## $ geom <POINT [m]> POINT (466401.1 4355852), POINT (46660…
## $ clump_id <fct> 471, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ clump_n_trees <int> 1, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 6…
## $ clump_n_trees_grp <fct> Individual, >25 trees, >25 trees, >25 trees,…
## $ clump_dist_m <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,…
## $ clump_n_trees_lb <fct> 1, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 2…
## $ clump_n_trees_ub <fct> 1, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, I…
## $ clump_nn_dist_m <dbl> NA, 4.031129, 4.756574, 3.905125, 4.596194, …
## $ clump_ostory_ht_m <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ clump_ostory_dbh_cm <dbl> 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12…let’s see a quick summary of the number of trees in each clump size
treetops_clustered %>%
sf::st_drop_geometry() %>%
dplyr::count(clump_n_trees_grp) %>%
dplyr::mutate(
pct_tot_trees = scales::percent(n/sum(n),accuracy=0.1)
) %>%
dplyr::rename(trees = n)## # A tibble: 6 × 3
## clump_n_trees_grp trees pct_tot_trees
## <fct> <int> <chr>
## 1 Individual 354 4.0%
## 2 2-4 trees 731 8.3%
## 3 5-9 trees 607 6.9%
## 4 10-15 trees 356 4.0%
## 5 16-25 trees 414 4.7%
## 6 >25 trees 6346 72.0%we can quickly look at the tree groups
plt_stand +
ggplot2::geom_sf(
data = treetops_clustered %>%
# rename the active geometry column to "geometry"
sf::st_set_geometry("geometry") %>%
dplyr::mutate(geometry = sf::st_buffer(geometry, dist = clump_dist_m/2))
, mapping = ggplot2::aes(fill = clump_id)
, color = NA
, show.legend = F
) +
ggplot2::scale_fill_manual(
values = c(
viridis::turbo( n = ceiling(length(unique(treetops_clustered$clump_id))/3) )
, viridis::viridis( n = ceiling(length(unique(treetops_clustered$clump_id))/3) )
, viridis::cividis( n = ceiling(length(unique(treetops_clustered$clump_id))/3) )
) %>% sample()
) 
Clump Polygons and Metrics
Churchill et al. (2016) provide instructions for implementing the clump identification (Plotkin et al. 2002) in ArcGIS:
Use the Buffer tool (in the Proximity toolset within the Analysis toolbox) to create a buffer of distance d/2, one half the inter-tree distance, around each point. This quantity d/2 is meant to approximate the crown radius of a “typical” overstory tree. Set the Dissolve Type option to ALL, which dissolves overlapping buffers, creating a reduced set of spatially non-overlapping polygons stored as a multipart polygon feature…Sanchez Meador et al. (2011) provide some useful examples of how clump attributes can be summarized…The method described here can be modified to use measured or modeled crown radii for each tree in place of d/2 (p.36)
# create function to pass a return from st_clump_points/get_tree_clumps and create clump polygons with summary stats
get_clump_summary <- function(get_tree_clumps_ans){
# get clump_dist_m
clump_dist_m <- min(get_tree_clumps_ans$clump_dist_m, na.rm = T)
# create clump polys and summary
clump_polys_temp <-
get_tree_clumps_ans %>%
dplyr::ungroup() %>%
# rename the active geometry column to "geometry"
sf::st_set_geometry("geometry") %>%
dplyr::mutate(geometry = sf::st_buffer(geometry, dist = clump_dist_m/2)) %>%
dplyr::group_by(clump_id, clump_n_trees_grp, clump_n_trees_lb, clump_n_trees_ub) %>%
dplyr::summarise(
# union buffered tree points
geometry = sf::st_union(geometry)
# summary metrics
, n_trees = dplyr::n_distinct(treeID)
, mean_dbh_cm = mean(dbh_cm, na.rm = T)
, min_dbh_cm = min(dbh_cm, na.rm = T)
, max_dbh_cm = max(dbh_cm, na.rm = T)
, mean_tree_height_m = mean(tree_height_m, na.rm = T)
, min_tree_height_m = min(tree_height_m, na.rm = T)
, max_tree_height_m = max(tree_height_m, na.rm = T)
, loreys_height_m = sum(basal_area_m2*tree_height_m, na.rm = T) / sum(basal_area_m2, na.rm = T)
, basal_area_m2 = sum(basal_area_m2, na.rm = T)
, sum_dbh_cm_sq = sum(dbh_cm^2, na.rm = T)
, .groups = "drop_last"
) %>%
dplyr::ungroup() %>%
sf::st_make_valid() %>%
dplyr::mutate(
clump_area_ha = sf::st_area(geometry) %>% as.numeric() %>% `/`(10000)
, trees_per_ha = (n_trees/clump_area_ha)
, basal_area_m2_per_ha = (basal_area_m2/clump_area_ha)
, pct_stand_basal_area = basal_area_m2/sum(basal_area_m2)
, pct_stand_n_trees = n_trees/sum(n_trees)
, qmd_cm = sqrt(sum_dbh_cm_sq/n_trees)
) %>%
dplyr::select(-c(sum_dbh_cm_sq)) %>%
# convert to imperial units
calc_imperial_units_fn() %>%
dplyr::mutate(clump_dist_m = clump_dist_m)
# calculate distance between clumps
clump_polys_temp <- clump_polys_temp %>%
dplyr::mutate(
nearest = sf::st_nearest_feature(clump_polys_temp)
, distance_nearest_clump_m = sf::st_distance(
clump_polys_temp
, clump_polys_temp[nearest,]
, by_element=TRUE
) %>%
as.numeric()
) %>%
dplyr::select(-c(nearest))
# return
return(clump_polys_temp)
}get the clump summary data using get_clump_summary()
# get it
clump_polys <- get_clump_summary(get_tree_clumps_ans = treetops_clustered)
# what?
clump_polys %>% dplyr::glimpse()## Rows: 824
## Columns: 34
## $ clump_id <fct> 471, 472, 473, 474, 475, 476, 477, 478, 479, …
## $ clump_n_trees_grp <fct> Individual, Individual, Individual, Individua…
## $ clump_n_trees_lb <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ clump_n_trees_ub <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ geometry <POLYGON [m]> POLYGON ((466404.1 4355852,..., POLYG…
## $ n_trees <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ mean_dbh_cm <dbl> 34.92759, 34.23414, 34.03887, 33.66477, 33.16…
## $ min_dbh_cm <dbl> 34.92759, 34.23414, 34.03887, 33.66477, 33.16…
## $ max_dbh_cm <dbl> 34.92759, 34.23414, 34.03887, 33.66477, 33.16…
## $ mean_tree_height_m <dbl> 24.723, 23.115, 22.135, 21.374, 21.033, 20.72…
## $ min_tree_height_m <dbl> 24.723, 23.115, 22.135, 21.374, 21.033, 20.72…
## $ max_tree_height_m <dbl> 24.723, 23.115, 22.135, 21.374, 21.033, 20.72…
## $ loreys_height_m <dbl> 24.723, 23.115, 22.135, 21.374, 21.033, 20.72…
## $ basal_area_m2 <dbl> 0.09581358, 0.09204683, 0.09099974, 0.0890104…
## $ clump_area_ha <dbl> 0.002826142, 0.002826142, 0.002826142, 0.0028…
## $ trees_per_ha <dbl> 353.8393, 353.8393, 353.8393, 353.8393, 353.8…
## $ basal_area_m2_per_ha <dbl> 33.90261, 32.56979, 32.19929, 31.49541, 30.56…
## $ pct_stand_basal_area <dbl> 0.0001779436, 0.0001709480, 0.0001690034, 0.0…
## $ pct_stand_n_trees <dbl> 0.0001135332, 0.0001135332, 0.0001135332, 0.0…
## $ qmd_cm <dbl> 34.92759, 34.23414, 34.03887, 33.66477, 33.16…
## $ mean_dbh_in <dbl> 13.76147, 13.48825, 13.41132, 13.26392, 13.06…
## $ min_dbh_in <dbl> 13.76147, 13.48825, 13.41132, 13.26392, 13.06…
## $ max_dbh_in <dbl> 13.76147, 13.48825, 13.41132, 13.26392, 13.06…
## $ qmd_in <dbl> 13.76147, 13.48825, 13.41132, 13.26392, 13.06…
## $ mean_tree_height_ft <dbl> 81.09144, 75.81720, 72.60280, 70.10672, 68.98…
## $ min_tree_height_ft <dbl> 81.09144, 75.81720, 72.60280, 70.10672, 68.98…
## $ max_tree_height_ft <dbl> 81.09144, 75.81720, 72.60280, 70.10672, 68.98…
## $ loreys_height_ft <dbl> 81.09144, 75.81720, 72.60280, 70.10672, 68.98…
## $ basal_area_ft2_per_ac <dbl> 147.7815, 141.9717, 140.3567, 137.2885, 133.2…
## $ trees_per_ac <dbl> 143.3049, 143.3049, 143.3049, 143.3049, 143.3…
## $ clump_area_ac <dbl> 0.006983396, 0.006983396, 0.006983396, 0.0069…
## $ basal_area_ft2 <dbl> 1.0313374, 0.9907921, 0.9795212, 0.9581088, 0…
## $ clump_dist_m <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, …
## $ distance_nearest_clump_m <dbl> 0.006544238, 0.267026450, 0.006544237, 1.6041…we should have the same number of polygons in the summary data as
identified in the get_tree_clumps() function
# do these numbers match
identical(
# clump
nrow(clump_polys)
# clumps in tree list data
, treetops_clustered %>% dplyr::distinct(clump_id) %>% nrow()
)## [1] TRUEplot the clump group size
plt_stand +
ggplot2::geom_sf(
data = clump_polys
, mapping = ggplot2::aes(fill = clump_n_trees_grp)
, color = NA
) +
harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
labs(fill = "") +
ggplot2::guides(size = "none", fill = ggplot2::guide_legend(override.aes = list(size = 5)))
ggplot2::ggsave("../data/02_tree_clumps.jpg", width = 7, height = 5)Clump Spacing
See Churchill et al. (2016) Figure 4 (p.10) and Matonis and Binkley (2018) who “calculated coverage of mosaic-meadows (percentage of stand > 6 m from overstory trees)” (p. 124)
Since we already buffered the tree points to approximate the crown radius, we’ll continue to use our \(d/2\) where \(d\) is maximum distance between trees for determining tree clumps and is meant to approximate the crown radius of a “typical” overstory tree
# create function to pass a return from get_clump_summary() and get a distance raster
get_clump_dist_rast <- function(get_clump_summary_ans, stand_sf){
# get clump_dist_m
clump_dist_m <- min(get_clump_summary_ans$clump_dist_m, na.rm = T)
# rasterize the clump polygons and then calculate distance between clumps as raster
dist_rast <-
terra::rasterize(
x = get_clump_summary_ans %>% terra::vect()
, y = get_clump_summary_ans %>%
terra::vect() %>%
terra::rast(res = 0.2)
) %>%
terra::distance() %>%
# crop it to stand extent
terra::crop(
stand_sf %>%
terra::vect()
, mask = T
)
######### part 2
# now create openings vector data
openings_vect <-
dist_rast %>%
terra::classify(rcl = c(clump_dist_m/2,Inf), others = NA, include.lowest = T) %>%
terra::as.polygons(na.rm = T) %>%
sf::st_as_sf() %>%
sf::st_cast("POLYGON") %>%
dplyr::mutate(layer = dplyr::row_number()) %>%
dplyr::mutate(
openining_area_m2 = sf::st_area(geometry) %>% as.numeric()
, clump_dist_m = clump_dist_m
)
# return
return(list(dist_rast = dist_rast, openings_vect = openings_vect))
}
# get it
get_clump_dist_rast_ans <- get_clump_dist_rast(clump_polys, stand_sf = stand_sf)what is in it?
names(get_clump_dist_rast_ans)## [1] "dist_rast" "openings_vect"just view the distance raster
terra::plot(get_clump_dist_rast_ans$dist_rast, axes = F, main = "Distance to Nearest Clump (m)")
plot the distance raster and openings vector data we just got with overlaid tree clumps and tree points
plt_fnl_1 <-
# ggplot2::ggplot() +
plt_stand +
# distance
ggplot2::geom_tile(
data = get_clump_dist_rast_ans$dist_rast %>% terra::aggregate(2, cores = 4) %>% terra::as.data.frame(xy = T) %>% dplyr::rename(f=3)
, mapping = ggplot2::aes(x=x, y=y, fill = f)
) +
harrypotter::scale_fill_hp(
"mischief"
, na.value = "transparent"
, name = "distance to\nnearest tree (m)"
, limits = c(0,22)
) +
# openings
ggplot2::geom_sf(
data = get_clump_dist_rast_ans$openings_vect
, mapping = ggplot2::aes(color = openining_area_m2), fill = NA
) +
ggplot2::scale_color_gradient(
low = "gray77", high = "gray11"
, labels = scales::comma_format(accuracy = 1)
, name = latex2exp::TeX("opening\narea ($\\m^2$)")
, limits = c(0,11100)
) +
# clumps
ggnewscale::new_scale_fill() +
ggplot2::geom_sf(
data = clump_polys
, mapping = ggplot2::aes(fill = clump_n_trees_grp)
, color = NA
) +
harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
# tree points
ggplot2::geom_sf(data = treetops_clustered, color = "gray88", shape = ".") +
ggplot2::theme_void() +
ggplot2::theme(
legend.position = "right"
, legend.direction = "vertical"
, plot.subtitle = ggplot2::element_text(size = 10, hjust = 0.5, face = "bold")
, legend.title= ggplot2::element_text(size=8)
, legend.text= ggplot2::element_text(size = 7)
)plot it
# plot
plt_fnl_1
# save it
ggplot2::ggsave("../data/03_clumps_dist.jpg", width = 7, height = 5)highlight the openings
plt_open_1 <-
plt_stand +
# clumps
ggplot2::geom_sf(
data = clump_polys
, mapping = ggplot2::aes(fill = clump_n_trees_grp)
, color = NA
) +
harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
# openings
ggnewscale::new_scale_fill() +
ggplot2::geom_sf(
data = get_clump_dist_rast_ans$openings_vect
, mapping = ggplot2::aes(fill = openining_area_m2), color = NA
) +
ggplot2::scale_fill_gradient(
low = "gray77", high = "gray11"
, labels = scales::comma_format(accuracy = 1)
, name = latex2exp::TeX("opening\narea ($\\m^2$)")
, limits = c(0,11100)
) +
# tree points
ggplot2::geom_sf(data = treetops_clustered, color = "gray88", shape = ".") +
ggplot2::theme_void() +
ggplot2::theme(
legend.position = "right"
, legend.direction = "vertical"
, plot.subtitle = ggplot2::element_text(size = 10, hjust = 0.5, face = "bold")
, legend.title= ggplot2::element_text(size=8)
, legend.text= ggplot2::element_text(size = 7)
)
# plot
plt_open_1
ggplot2::ggsave("../data/04_clumps_openings.jpg", width = 7, height = 5)combine them with patchwork
plt_fnl_1 + (plt_open_1 + theme(legend.position = "none")) 
# &
# ggplot2::theme(
# plot.subtitle = ggplot2::element_text(size = 10, hjust = 0.5, face = "bold")
# , legend.title=ggplot2::element_text(size=8)
# , legend.text=ggplot2::element_text(size = 7)
# )
ggplot2::ggsave("../data/05_clumps_combined.jpg", width = 10, height = 5.5)Clump Metrics
create a function to summarize by the tree clump grouping variable (e.g. “Individual”,“2-4 trees”,“5-9 trees”,“10-15 trees”,…)
# create a function to summarize by number of tree clump grouping
get_clump_n_trees_grp_summary = function(get_tree_clumps_ans, get_clump_summary_ans, stand_sf){
# get area of harvest unit
#...will use this area in the area calculations such that...
#...TPA = trees in a certain group size across the whole stand area
harvest_area_ha <- stand_sf %>%
sf::st_union() %>%
sf::st_as_sf() %>%
dplyr::mutate(harvest_area_ha = sf::st_area(.) %>% as.numeric() %>% `/`(10000)) %>%
dplyr::pull(harvest_area_ha) %>%
.[1]
# collapse and calculate silv metrics
dta <- get_tree_clumps_ans %>%
sf::st_drop_geometry() %>%
dplyr::mutate(stand_area_ha = harvest_area_ha) %>%
dplyr::group_by(stand_area_ha,clump_n_trees_grp, clump_n_trees_lb, clump_n_trees_ub) %>%
dplyr::summarise(
# summary metrics
n_trees = dplyr::n_distinct(treeID)
, mean_dbh_cm = mean(dbh_cm, na.rm = T)
, mean_tree_height_m = mean(tree_height_m, na.rm = T)
, loreys_height_m = sum(basal_area_m2*tree_height_m, na.rm = T) / sum(basal_area_m2, na.rm = T)
, basal_area_m2 = sum(basal_area_m2, na.rm = T)
, sum_dbh_cm_sq = sum(dbh_cm^2, na.rm = T)
, .groups = "drop_last"
) %>%
dplyr::ungroup() %>%
# attach clump area
dplyr::left_join(
get_clump_summary_ans %>%
sf::st_drop_geometry() %>%
dplyr::group_by(clump_n_trees_grp) %>%
dplyr::summarise(
clump_area_ha = sum(clump_area_ha)
, stand_n_clumps = dplyr::n()
, .groups = "drop_last"
) %>%
dplyr::ungroup()
, by = dplyr::join_by(clump_n_trees_grp)
) %>%
dplyr::mutate(
trees_per_ha = (n_trees/stand_area_ha) # (n_trees/clump_area_ha) ... this was not right
, basal_area_m2_per_ha = (basal_area_m2/stand_area_ha) # (basal_area_m2/clump_area_ha) ... this was not right
, qmd_cm = sqrt(sum_dbh_cm_sq/n_trees)
# stand calcs
, stand_trees_per_ha = sum(n_trees)/stand_area_ha
, stand_basal_area_m2 = sum(basal_area_m2)
, stand_basal_area_m2_per_ha = sum(basal_area_m2)/stand_area_ha
, pct_stand_basal_area = basal_area_m2/stand_basal_area_m2
, pct_stand_n_trees = n_trees/sum(n_trees)
, stand_qmd_cm = sqrt(sum(get_tree_clumps_ans$dbh_cm^2, na.rm = T)/sum(n_trees))
) %>%
dplyr::select(-c(sum_dbh_cm_sq)) %>%
# convert to imperial units
calc_imperial_units_fn()
# return
return(dta)
}call it
# call it
clump_n_trees_grp_summary <- get_clump_n_trees_grp_summary(
get_tree_clumps_ans = treetops_clustered
, get_clump_summary_ans = clump_polys
, stand_sf = stand_sf
)
# what?
clump_n_trees_grp_summary %>% dplyr::glimpse()## Rows: 6
## Columns: 33
## $ stand_area_ha <dbl> 30.4578, 30.4578, 30.4578, 30.4578, 30.457…
## $ clump_n_trees_grp <fct> Individual, 2-4 trees, 5-9 trees, 10-15 tr…
## $ clump_n_trees_lb <fct> 1, 2, 5, 10, 16, 26
## $ clump_n_trees_ub <fct> 1, 4, 9, 15, 25, Inf
## $ n_trees <int> 354, 731, 607, 356, 414, 6346
## $ mean_dbh_cm <dbl> 27.69943, 27.64553, 26.96353, 26.74255, 26…
## $ mean_tree_height_m <dbl> 15.36712, 15.19598, 14.60823, 14.55307, 14…
## $ loreys_height_m <dbl> 16.79629, 16.48573, 15.94476, 16.09756, 15…
## $ basal_area_m2 <dbl> 22.22600, 45.55548, 36.09377, 21.00662, 24…
## $ clump_area_ha <dbl> 1.0004541, 1.7604321, 1.3383869, 0.7875838…
## $ stand_n_clumps <int> 354, 276, 98, 30, 20, 46
## $ trees_per_ha <dbl> 11.62264, 24.00042, 19.92921, 11.68830, 13…
## $ basal_area_m2_per_ha <dbl> 0.7297307, 1.4956915, 1.1850417, 0.6896957…
## $ qmd_cm <dbl> 28.27379, 28.16869, 27.51546, 27.40994, 27…
## $ stand_trees_per_ha <dbl> 289.187, 289.187, 289.187, 289.187, 289.18…
## $ stand_basal_area_m2 <dbl> 538.4492, 538.4492, 538.4492, 538.4492, 53…
## $ stand_basal_area_m2_per_ha <dbl> 17.67853, 17.67853, 17.67853, 17.67853, 17…
## $ pct_stand_basal_area <dbl> 0.04127779, 0.08460496, 0.06703281, 0.0390…
## $ pct_stand_n_trees <dbl> 0.04019074, 0.08299273, 0.06891462, 0.0404…
## $ stand_qmd_cm <dbl> 27.89901, 27.89901, 27.89901, 27.89901, 27…
## $ mean_dbh_in <dbl> 10.91357, 10.89234, 10.62363, 10.53656, 10…
## $ qmd_in <dbl> 11.13987, 11.09846, 10.84109, 10.79952, 10…
## $ stand_qmd_in <dbl> 10.99221, 10.99221, 10.99221, 10.99221, 10…
## $ mean_tree_height_ft <dbl> 50.40415, 49.84282, 47.91499, 47.73406, 46…
## $ loreys_height_ft <dbl> 55.09184, 54.07318, 52.29883, 52.80000, 50…
## $ basal_area_ft2_per_ac <dbl> 3.180896, 6.519719, 5.165597, 3.006383, 3.…
## $ stand_basal_area_ft2_per_ac <dbl> 77.06072, 77.06072, 77.06072, 77.06072, 77…
## $ trees_per_ac <dbl> 4.707168, 9.720169, 8.071330, 4.733762, 5.…
## $ stand_trees_per_ac <dbl> 117.1207, 117.1207, 117.1207, 117.1207, 11…
## $ stand_area_ac <dbl> 75.26123, 75.26123, 75.26123, 75.26123, 75…
## $ clump_area_ac <dbl> 2.472122, 4.350028, 3.307154, 1.946119, 2.…
## $ basal_area_ft2 <dbl> 239.2406, 490.3592, 388.5133, 226.1152, 25…
## $ stand_basal_area_ft2 <dbl> 5795.868, 5795.868, 5795.868, 5795.868, 57…summary table
# table it
clump_n_trees_grp_summary %>%
dplyr::select(
clump_n_trees_grp, n_trees
, mean_dbh_in
, qmd_in
, mean_tree_height_ft
, loreys_height_ft
, trees_per_ac
, basal_area_ft2_per_ac, pct_stand_basal_area, pct_stand_n_trees
) %>%
dplyr::mutate(
dplyr::across(
.cols = c(pct_stand_basal_area, pct_stand_n_trees)
, .fns = ~ scales::percent(.x, accuracy = 1)
)
) %>%
kableExtra::kbl(
digits = 1
, escape = F
, caption = paste0("Overstory tree clump summary")
, col.names = c(
"", "trees"
, "mean<br>DBH (in)"
, "QMD (in)"
, "mean<br>Ht. (ft)"
, "Loreys<br>Ht. (ft)"
, "TPA"
, "BA<br>ft<sup>2</sup> ac<sup>-1</sup>"
, "Pct<br>stand BA"
, "Pct<br>stand trees"
)
) %>%
kableExtra::kable_styling()| trees |
mean DBH (in) |
QMD (in) |
mean Ht. (ft) |
Loreys Ht. (ft) |
TPA |
BA ft2 ac-1 |
Pct stand BA |
Pct stand trees |
|
|---|---|---|---|---|---|---|---|---|---|
| Individual | 354 | 10.9 | 11.1 | 50.4 | 55.1 | 4.7 | 3.2 | 4% | 4% |
| 2-4 trees | 731 | 10.9 | 11.1 | 49.8 | 54.1 | 9.7 | 6.5 | 8% | 8% |
| 5-9 trees | 607 | 10.6 | 10.8 | 47.9 | 52.3 | 8.1 | 5.2 | 7% | 7% |
| 10-15 trees | 356 | 10.5 | 10.8 | 47.7 | 52.8 | 4.7 | 3.0 | 4% | 4% |
| 16-25 trees | 414 | 10.5 | 10.7 | 46.9 | 50.8 | 5.5 | 3.5 | 4% | 5% |
| >25 trees | 6346 | 10.9 | 11.0 | 48.6 | 51.4 | 84.4 | 55.7 | 72% | 72% |
ICO Implementation
Churchill et al. (2016) describe the full process for implementing the ICO approach in The ICO Approach to Quantifying and Restoring Forest Spatial Pattern: Implementation Guide in which the authors lay out the prescription development process:
- Identify skips and other special treatment areas
- Consider the need for openings
- Determine the stand average density target
- Determine the appropriate distance to define clumps
- Obtain targets for clump proportions
- Select target clump proportions for your stand
- Generate clump targets for the whole unit
- Combine clump and opening targets with leave tree criteria into marking guidelines
The objective here is to: 1) provide the manager with the current conditions (completed above); 2) take the “targets” as set by the manager (steps 3, 5, 6, 7); 3) create the prescription with the leave tree marking.
Let’s implement this prescription development process with our UAS tree list
3. Determine the stand average density target
Step 3 in Churchill et al. (2016):
An average BA, TPA, or SDI target for the stand should be selected that is appropriate for the species, structure, site conditions, and management objectives. Expected mortality from prescribed fire should be factored in. Stand average targets can come from historical reference stands, plant association based stocking guides, density management tools, or a combination of both (see Franklin et al. (2013) for a full discussion of setting density targets). In dry forests, the number and size of old trees must be accounted in setting the density target. To use the ICO method, the target must be converted to TPA (see Table 1). A lower diameter cutoff also needs to be specified for the TPA target. This should be the lower limit in the contract or cutting guidelines given to the marking crew or contractor. (p.11)
this is what Table 1 looks like with TPA values are derived from the formula:
\[ TPA = \frac{BA}{QMD^{2} \times 0.005454} \]
# function to get tpa from ba and qmd
get_tpa <- function(ba_ft2_ac, qmd_in){
tpa <- round(ba_ft2_ac/((qmd_in^2)*0.005454))
return(tpa)
}
# table it
tidyr::crossing(
ba = seq(40,200,20)
, qmd = seq(8,20,2)
) %>%
dplyr::mutate(
tpa = get_tpa(ba,qmd)
) %>%
tidyr::pivot_wider(names_from = ba, values_from = tpa) %>%
dplyr::mutate(l = "QMD (in)") %>%
dplyr::relocate(l) %>%
kableExtra::kbl(
col.names = c(".","", seq(40,200,20))
, escape = F
, caption = "Basal Area and QMD to TPA conversion chart"
) %>%
kableExtra::add_header_above(
c("","", "Basal Area (ft<sup>2</sup> ac<sup>-1</sup>)"=length(seq(40,200,20)))
, escape = F
) %>%
kableExtra::kable_styling() %>%
kableExtra::column_spec(1:2, bold = T) %>%
kableExtra::collapse_rows(columns = 1, valign = "middle")| . | 40 | 60 | 80 | 100 | 120 | 140 | 160 | 180 | 200 | |
|---|---|---|---|---|---|---|---|---|---|---|
| QMD (in) | 8 | 115 | 172 | 229 | 286 | 344 | 401 | 458 | 516 | 573 |
| 10 | 73 | 110 | 147 | 183 | 220 | 257 | 293 | 330 | 367 | |
| 12 | 51 | 76 | 102 | 127 | 153 | 178 | 204 | 229 | 255 | |
| 14 | 37 | 56 | 75 | 94 | 112 | 131 | 150 | 168 | 187 | |
| 16 | 29 | 43 | 57 | 72 | 86 | 100 | 115 | 129 | 143 | |
| 18 | 23 | 34 | 45 | 57 | 68 | 79 | 91 | 102 | 113 | |
| 20 | 18 | 28 | 37 | 46 | 55 | 64 | 73 | 83 | 92 |
Current Stand Conditions
For determining targets, the silviculturist needs to know the current conditions. Provide the current stand conditions based on the UAS tree list for the selected stand that are required to set the targets:
- Current BA
- Current QMD
- Current proportion of trees by clump size
clump_n_trees_grp_summary %>%
dplyr::select(clump_n_trees_grp, pct_stand_n_trees) %>%
dplyr::mutate(
pct_stand_n_trees = scales::percent(pct_stand_n_trees,accuracy = 1)
) %>%
tidyr::pivot_wider(names_from = clump_n_trees_grp, values_from = pct_stand_n_trees) %>%
kableExtra::kable(
caption = paste0(
"Current stand BA (ft<sup>2</sup> ac<sup>-1</sup>): "
, clump_n_trees_grp_summary$stand_basal_area_ft2_per_ac[1] %>% scales::number(accuracy = 0.1)
, "<br>Current stand QMD (in): "
, clump_n_trees_grp_summary$stand_qmd_in[1] %>% scales::number(accuracy = 0.1)
, "<br>Current stand TPA: "
, clump_n_trees_grp_summary$stand_trees_per_ac[1] %>% scales::number(accuracy = 1)
)
, escape = F
, digits = 1
) %>%
kableExtra::kable_styling() %>%
kableExtra::footnote(general = "values are the percent of trees in each clump size")| Individual | 2-4 trees | 5-9 trees | 10-15 trees | 16-25 trees | >25 trees |
|---|---|---|---|---|---|
| 4% | 8% | 7% | 4% | 5% | 72% |
| Note: | |||||
| values are the percent of trees in each clump size |
5. Obtain targets for clump proportions
Step 5 in Churchill et al. (2016):
ICO prescriptions are based on a target proportion of trees in different sized clumps within a stand. Proportions are just the percentage of trees, or TPA, that are in different sized clumps. Basal area proportions can be used, but we have found TPA targets to be more straightforward to use. Ideally, a table summarizing clump proportions for a range of reference conditions in your area is available (Table 2). If not, instructions for developing one are provided in section VI. (p.12)
Section VI of Churchill et al. (2016) notes that
reference spatial information may already be available and summarized in a way that it can be directly incorporated into ICO prescriptions. Such data exist and have been published for areas in Arizona (Abella and Denton 2009, Sánchez Meador et al. 2011), the eastern Washington Cascades (Churchill et al. 2013), the northern Rockies (Larson et al. 2012), and the Sierra Nevada (Lydersen et al. 2013). Reference datasets for using ICO in other forest types, such as coastal Douglas-fir or Pacific silver fir, also exist (Larson and Churchill 2008). (p.28)
Table 2 is:
6. Select target clump proportions for your stand
Now set the desired BA, QMD, and proportion of trees in each clump size:
########################################################################################
########################################################################################
# desired BA, QMD, and proportion of trees in each clump size
########################################################################################
########################################################################################
# desired BA
# target_ba_ft2_per_ac <- 9 # cannot be > current BA
target_ba_ft2_per_ac <- 50 # cannot be > current BA
# desired QMD
target_qmd_in <- 12
# desired proportion (%) of trees in each clump size
# !cannot be create larger proportion of ">25 trees" clump as this would require adding trees...
# c("Individual", "2-4 trees", "5-9 trees", "10-15 trees","16-25 trees",">25 trees")
# c(.18, .33, .24, .10, .15)
# target_pcts <- c(66, 33, 1, 0, 0, 0)
target_pcts <- c(25, 40, 25, 10, 0, 0)
########################################################################################
########################################################################################
# desired BA, QMD, and proportion of trees in each clump size
########################################################################################
########################################################################################Check set up and define data with targets based on Churchill et al. (2016) and their prescription worksheet Excel file to help develop ICO prescriptions
get_target_check_prescription = function(
clump_n_trees_grp_summary_dta
, target_ba_ft2_per_ac = NA
, target_qmd_in = NA
, target_pcts = NA
){
if(
is.na(target_ba_ft2_per_ac) | is.na(target_qmd_in) | any(is.na(target_pcts))
){
stop("must set all of the function parameters:\n`target_ba_ft2_per_ac`, `target_qmd_in`, and `target_pcts`")
}
#############################################
# check target BA and TPA
#############################################
if(as.numeric(target_ba_ft2_per_ac)>clump_n_trees_grp_summary_dta$stand_basal_area_ft2_per_ac[1]){
stop(
"target BA in `target_ba_ft2_per_ac` of "
, round(as.numeric(target_ba_ft2_per_ac),1), " is greater than current BA of "
, clump_n_trees_grp_summary_dta$stand_basal_area_ft2_per_ac[1] %>% round(1)
)
}
if(
get_tpa(target_ba_ft2_per_ac, target_qmd_in)>clump_n_trees_grp_summary_dta$stand_trees_per_ac[1]
){
stop(
"target TPA based on `target_ba_ft2_per_ac` and `target_qmd_in` of "
, round(as.numeric(get_tpa(target_ba_ft2_per_ac, target_qmd_in)),1), " is greater than current TPA of "
, clump_n_trees_grp_summary_dta$stand_trees_per_ac[1] %>% round(1)
, "\n adjust `target_ba_ft2_per_ac` and/or `target_qmd_in` to get valid TPA"
)
}
# if(sum(target_pcts)!=1){stop("`target_pcts` must sum to 1 (100%)")}
#############################################
# define data with current and target
# ... this is "smart" in that percentages are adj based on:
# ... 0) are there missing targets?
# ... ... if < 6 numbers provided in `target_pcts` then the largest tree groups get targets of 0
# ... 1) do targets sum to 1?
# ... ... if not trees are distributed proportionally based on targets provided and trees available
# ... 2) is target in largest clump size > current conditions?
# ... ... if yes, target is set to current condition
# ... 3) is target listed in clump size > current largest clump with trees?
# ... ... if yes, target for largest clump size is shifted to current largest clump with trees
#############################################
# define clump group data even if levels are missing from data
clump_n_trees_grp_df <-
dplyr::tibble(
# get all possible labels from the factor levels
clump_n_trees_grp = levels(clump_n_trees_grp_summary_dta$clump_n_trees_grp)
, clump_n_trees_lb = levels(clump_n_trees_grp_summary_dta$clump_n_trees_lb)
, clump_n_trees_ub = levels(clump_n_trees_grp_summary_dta$clump_n_trees_ub)
) %>%
# ensure they are still factors
dplyr::mutate(
dplyr::across(
dplyr::everything()
, ~factor(.x, levels = .x)
)
# get numeric values
, min_clump_n_trees = clump_n_trees_lb %>% as.character() %>% as.numeric()
, max_clump_n_trees = clump_n_trees_ub %>% as.character() %>% as.numeric()
) %>%
# calculate the median to use as the "average" of the clump size
dplyr::rowwise() %>%
dplyr::mutate(
mean_clump_n_trees = median(c( # note the c()
min_clump_n_trees
, ifelse( is.infinite(max_clump_n_trees), ceiling(min_clump_n_trees*1.3), max_clump_n_trees )
)) %>%
floor()
# c(1,3,7,12,20,30)
) %>%
dplyr::ungroup()
target_data <-
# create data for joining if missing clump groups
clump_n_trees_grp_df %>%
dplyr::mutate(
stand_area_ac = clump_n_trees_grp_summary_dta$stand_area_ac[1]
) %>%
dplyr::left_join(
clump_n_trees_grp_summary_dta %>%
dplyr::ungroup() %>%
dplyr::select(clump_n_trees_grp, pct_stand_n_trees, stand_n_clumps)
, by = "clump_n_trees_grp"
) %>%
dplyr::mutate(
pct_stand_n_trees = dplyr::coalesce(pct_stand_n_trees,0)
, stand_n_clumps = dplyr::coalesce(stand_n_clumps,0)
) %>%
# add targets
dplyr::bind_cols(
pct_stand_n_trees_target = c(
as.numeric(target_pcts)
, rep(0,times=nrow(clump_n_trees_grp_df))
)[1:nrow(clump_n_trees_grp_df)] # pad target with 0's
) %>%
# adjust target based on difference from 1
dplyr::mutate(
pct_stand_n_trees_target = pct_stand_n_trees_target*(1/sum(pct_stand_n_trees_target))
# largest clump size with trees
, largest_w_trees = max(
ifelse(
dplyr::coalesce(pct_stand_n_trees)>0
,min_clump_n_trees
,NA
)
,na.rm = T
)
, largest_w_trees_target = max(
ifelse(
dplyr::coalesce(pct_stand_n_trees_target)>0
,min_clump_n_trees
,NA
)
,na.rm = T
)
) %>%
# move target for largest clump size to the largest current clump size
dplyr::mutate(
pct_stand_n_trees_target = dplyr::case_when(
min_clump_n_trees==largest_w_trees &
largest_w_trees_target>largest_w_trees ~ max(
ifelse(min_clump_n_trees==largest_w_trees_target,pct_stand_n_trees_target,0)
)
, T ~ pct_stand_n_trees_target
)
) %>%
# adjust target based on current conditions
dplyr::mutate(
pct_stand_n_trees_target = dplyr::case_when(
min_clump_n_trees>largest_w_trees &
pct_stand_n_trees_target > 0 ~ 0
, min_clump_n_trees==largest_w_trees &
pct_stand_n_trees_target > pct_stand_n_trees ~ pct_stand_n_trees
, T ~ pct_stand_n_trees_target
)
) %>%
# finally, re-scale again based on adjustments
dplyr::mutate(
pct_stand_n_trees_target = dplyr::case_when(
min_clump_n_trees==largest_w_trees ~ pct_stand_n_trees_target
, T ~ pct_stand_n_trees_target * (
# pct remaining to scale to
(1-max(ifelse(min_clump_n_trees==largest_w_trees,pct_stand_n_trees_target,0))) /
# current pct remaining total allocated
sum(
ifelse(min_clump_n_trees!=largest_w_trees,pct_stand_n_trees_target,0))
)
)
) %>%
# add other targets
dplyr::rename(pct_stand_n_trees_current = pct_stand_n_trees) %>%
dplyr::mutate(
stand_trees_per_ac_current = clump_n_trees_grp_summary_dta$stand_trees_per_ac[1]
, stand_trees_per_ac_target = dplyr::coalesce(get_tpa(target_ba_ft2_per_ac, target_qmd_in),0)
, trees_per_acre_current = stand_trees_per_ac_current*pct_stand_n_trees_current
, trees_per_acre_target = stand_trees_per_ac_target*pct_stand_n_trees_target
, clumps_per_acre_current = trees_per_acre_current/mean_clump_n_trees
, clumps_per_acre_target = trees_per_acre_target/mean_clump_n_trees
, stand_n_clumps_current = stand_n_clumps
, stand_n_clumps_target = (clumps_per_acre_target*stand_area_ac) %>% round(0)
) %>%
dplyr::select(-c(tidyselect::starts_with("largest_w_trees"), stand_n_clumps))
# ????
# target_data %>% glimpse()
# issue warning about targets
if(
!all(
target_data$pct_stand_n_trees_target == c(
as.numeric(target_pcts)
, rep(0,times=nrow(clump_n_trees_grp_df))
)[1:nrow(clump_n_trees_grp_df)]
)
){
warning(
"proportion of trees in each clump size target `target_pcts` adjusted!!!"
, "\nfrom : ", paste(round(target_pcts,2),collapse = ",")
, "\nto : ", paste(round(target_data$pct_stand_n_trees_target,2),collapse = ",")
)
}
# return
return(target_data)
}
# call it
prescription_target_data <- get_target_check_prescription(
clump_n_trees_grp_summary
, target_ba_ft2_per_ac = target_ba_ft2_per_ac
, target_qmd_in = target_qmd_in
, target_pcts = target_pcts
)## Warning in get_target_check_prescription(clump_n_trees_grp_summary, target_ba_ft2_per_ac = target_ba_ft2_per_ac, : proportion of trees in each clump size target `target_pcts` adjusted!!!
## from : 25,40,25,10,0,0
## to : 0.25,0.4,0.25,0.1,0,0# what?
prescription_target_data %>% dplyr::glimpse()## Rows: 6
## Columns: 17
## $ clump_n_trees_grp <fct> Individual, 2-4 trees, 5-9 trees, 10-15 tre…
## $ clump_n_trees_lb <fct> 1, 2, 5, 10, 16, 26
## $ clump_n_trees_ub <fct> 1, 4, 9, 15, 25, Inf
## $ min_clump_n_trees <dbl> 1, 2, 5, 10, 16, 26
## $ max_clump_n_trees <dbl> 1, 4, 9, 15, 25, Inf
## $ mean_clump_n_trees <dbl> 1, 3, 7, 12, 20, 30
## $ stand_area_ac <dbl> 75.26123, 75.26123, 75.26123, 75.26123, 75.…
## $ pct_stand_n_trees_current <dbl> 0.04019074, 0.08299273, 0.06891462, 0.04041…
## $ pct_stand_n_trees_target <dbl> 0.25, 0.40, 0.25, 0.10, 0.00, 0.00
## $ stand_trees_per_ac_current <dbl> 117.1207, 117.1207, 117.1207, 117.1207, 117…
## $ stand_trees_per_ac_target <dbl> 64, 64, 64, 64, 64, 64
## $ trees_per_acre_current <dbl> 4.707168, 9.720169, 8.071330, 4.733762, 5.5…
## $ trees_per_acre_target <dbl> 16.0, 25.6, 16.0, 6.4, 0.0, 0.0
## $ clumps_per_acre_current <dbl> 4.7071680, 3.2400563, 1.1530472, 0.3944802,…
## $ clumps_per_acre_target <dbl> 16.0000000, 8.5333333, 2.2857143, 0.5333333…
## $ stand_n_clumps_current <dbl> 354, 276, 98, 30, 20, 46
## $ stand_n_clumps_target <dbl> 1204, 642, 172, 40, 0, 0# !!!!!!!!!!!!!!!!!!!! need to ensure sufficient trees in larger clump groups if increasing target for lower groups
### e.g. see the third record:
# $ pct_stand_n_trees_current <dbl> 0.28, 0.43, 0.19, 0.07, 0.00, 0.01
# $ pct_stand_n_trees_target <dbl> 0.40, 0.30, 0.28, 0.02, 0.00, 0.00current vs target
prescription_target_data %>%
dplyr::select(
clump_n_trees_grp
, tidyselect::starts_with("pct_stand_n_trees")
, tidyselect::starts_with("trees_per_acre_")
) %>%
dplyr::mutate(
dplyr::across(
tidyselect::starts_with("pct_stand_n_trees")
, ~ scales::percent(.x,accuracy = 1)
)
) %>%
kableExtra::kable(
caption = paste0(
"Current stand (ft<sup>2</sup> ac<sup>-1</sup>): "
, clump_n_trees_grp_summary$stand_basal_area_ft2_per_ac[1] %>% scales::number(accuracy = 0.1)
, "<br>Current stand QMD (in): "
, clump_n_trees_grp_summary$stand_qmd_in[1] %>% scales::number(accuracy = 0.1)
, "<br>Current stand TPA: "
, clump_n_trees_grp_summary$stand_trees_per_ac[1] %>% scales::number(accuracy = 1)
)
, escape = F
, digits = 1
, col.names = c(
"", rep(c("current","target"),2)
)
) %>%
kableExtra::kable_styling() %>%
kableExtra::add_header_above(
c(" "=1,"% Trees"=2, "TPA"=2)
)| current | target | current | target | |
|---|---|---|---|---|
| Individual | 4% | 25% | 4.7 | 16.0 |
| 2-4 trees | 8% | 40% | 9.7 | 25.6 |
| 5-9 trees | 7% | 25% | 8.1 | 16.0 |
| 10-15 trees | 4% | 10% | 4.7 | 6.4 |
| 16-25 trees | 5% | 0% | 5.5 | 0.0 |
| >25 trees | 72% | 0% | 84.4 | 0.0 |
8. Combine clump and opening targets with leave tree criteria into marking guidelines
Use our UAS tree list to generate the prescription:
- start with the largest clump size currently with trees
- cut trees to the next largest clump size until desired # clumps is reached
- repeat with each successive clump size through to individual tree selection
- if possible, cut in same clump until desired proportions are reached to minimize machine time
Function to cut clumps
put the entire process outlined immediately above into a function to use with the tree list to cut each clump and then we’ll select the desired proportion of clumps
define intermediate functions
define intermediate functions
##############################################
# working with sf LINESTRINGS
##############################################
# first two functions borrowed from https://github.com/metafor-ulaval/ALSroads/blob/main/R/line_tools.R
########
# Get heading of both ends of a line
########
st_ends_heading <- function(line){
M <- sf::st_coordinates(line)
i <- c(2, nrow(M) - 1)
j <- c(1, -1)
headings <- mapply(i, j, FUN = function(i, j) {
Ax = M[i-j,1]
Ay = M[i-j,2]
Bx = M[i,1]
By = M[i,2]
atan2(Ay-By, Ax-Bx)*180/pi
})
names(headings) <- c("head", "tail")
return(headings)
}
########
# extend the line on both ends
########
st_extend_line <- function(line, distance, end = "BOTH"){
if (!(end %in% c("BOTH", "HEAD", "TAIL")) | length(end) != 1) stop("'end' must be 'BOTH', 'HEAD' or 'TAIL'")
M <- sf::st_coordinates(line)[,-3]
keep <- !(end == c("TAIL", "HEAD"))
ends <- c(1, nrow(M))[keep]
headings <- st_ends_heading(line)[keep] / 180 * pi
distances <- if (length(distance) == 1) rep(distance, 2) else distance[1:2]
M[ends, 1:2] <- M[ends, 1:2] + distances[keep] * c(cos(headings), sin(headings))
newline <- sf::st_linestring(M)
# If input is sfc_LINESTRING and not sfg_LINESTRING
if (is.list(line)) newline <- sf::st_sfc(newline, crs = sf::st_crs(line))
return(newline)
}
########
# pass an sf dataframe of points and return a line between the farthest points
########
st_points_to_line <- function(pts, line_ext=0) {
# find farthest distance between points
dist_temp <- sf::st_distance(pts)
# get the points
f_pts_temp <-
pts %>%
dplyr::ungroup() %>%
dplyr::slice(
# get the farthest points from distance matrix
which(dist_temp == max(dist_temp), arr.ind = TRUE)[1,]
)
# draw a line between the farthest two points
f_line_temp <- f_pts_temp %>%
# convert to linestring
dplyr::ungroup() %>%
dplyr::summarise(n=dplyr::n()) %>%
sf::st_cast("LINESTRING") %>%
dplyr::ungroup() %>%
dplyr::select(-c(n))
# and apply the line extension
farthest_line <- st_extend_line(f_line_temp, distance = line_ext)
# return
return(farthest_line)
}
# st_points_to_line(ttops_temp %>% dplyr::slice_head(prop = .1), line_ext = 6) %>%
# ggplot() + geom_sf() + geom_sf(data = ttops_temp %>% dplyr::slice_head(prop = .1)) + theme_void()
########
# Function to calculate Euclidean distance between 2 points
########
st_euclidean_distance <- function(p1,p2) {
return(sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2))
}
########
# return a line perpendicular to current line
########
### https://stackoverflow.com/questions/56771058/perpendicular-lines-at-regular-intervals-along-lines-with-multiple-nodes
# Function to calculate 2 points on a line perpendicular to another defined by 2 points p1,p2
# For point at interval, which can be a proportion of the segment length, or a constant
st_perp_line <- function(interval=0.5, my_line, proportion=TRUE) {
# get end points of line
p1 <- my_line %>% sf::st_cast("POINT") %>% sf::st_coordinates() %>% .[1,]
p2 <- my_line %>% sf::st_cast("POINT") %>% sf::st_coordinates() %>% .[2,]
# get length of line to return equal length line
line_len <- sf::st_length(my_line) %>% as.numeric() %>% `/`(2)
# get crs of line
my_crs <- sf::st_crs(my_line)
# Calculate x and y distances
x_len <- p2[1] - p1[1]
y_len <- p2[2] - p1[2]
# If proportion calculate reference point from tot_length
if (proportion) {
point <- c(p1[1]+x_len*interval,p1[2]+y_len*interval)
}
# Else use the constant value
else {
tot_len <- st_euclidean_distance(p1,p2)
point <- c(p1[1]+x_len/tot_len*interval,p1[2]+y_len/tot_len*interval)
}
# Calculate the x and y distances from reference point to point on line line_len distance away
ref_len <- st_euclidean_distance(point,p2)
xn_len <- (line_len / ref_len) * (p2[1] - point[1])
yn_len <- (line_len / ref_len) * (p2[2] - point[2])
# fix for identical
if(identical(point,p2) && x_len>y_len){ # this works for horizontal line
xn_len <- line_len/2
yn_len <- 0
}else if(identical(point,p2) && x_len<y_len){ # this works for vertical line
xn_len <- 0
yn_len <- line_len/2
}
# Invert the x and y lengths and add/subtract from the refrence point
# ref_points <- rbind(point,c(point[1] + yn_len,point[2] - xn_len),c(point[1] - yn_len,point[2] + xn_len))
ref_points <- rbind(c(point[1] + yn_len,point[2] - xn_len),c(point[1] - yn_len,point[2] + xn_len))
# use the reference points to return a line
return(
ref_points %>%
dplyr::as_tibble() %>%
dplyr::rename_with(tolower) %>%
sf::st_as_sf(coords = c("x", "y"), crs = my_crs, remove = F) %>%
dplyr::summarise(n=dplyr::n()) %>%
sf::st_cast("LINESTRING") %>%
dplyr::ungroup() %>%
dplyr::select(-c(n))
)
}function to pass a clump
generate data using clumping functions above, pass that data, return tree list with cut/keep flag based on the target clump group size
# treetops_clustered %>%
# dplyr::filter( as.numeric(as.character(clump_n_trees_lb)) == prescription_target_data$min_clump_n_trees[3]) %>%
# dplyr::arrange(desc(clump_n_trees)) %>%
# # dplyr::slice(1) %>%
# # dplyr::pull(clump_id)
# dplyr::filter(clump_id == 2) %>%
# dplyr::select(treeID, tidyselect::starts_with("clump_")) %>%
# dplyr::glimpse()
# levels(treetops_clustered$clump_n_trees_lb) %>% as.character() %>% as.numeric()
# c(
# 0 # always pad with 0
# , levels(treetops_clustered$clump_n_trees_ub) %>%
# as.character() %>%
# as.numeric() %>%
# unique() %>%
# purrr::discard(~ is.null(.x) || is.na(.x))
# , Inf # always pad with Inf
# ) %>%
# unique() %>%
# sort()
# function to cut a single clump
# as.numeric(as.character(levels(treetops_clustered$clump_n_trees_lb)))
# let's define a function to work with the target and the factors from st_clump_points
check_clump_breaks <- function(st_clump_points_ans) {
# check cols
cloud2trees:::check_df_cols_all_missing(
st_clump_points_ans
, col_names = c("clump_n_trees_grp","clump_n_trees_lb", "clump_n_trees_ub","clump_dist_m")
, all_numeric = F, check_vals_missing = F
)
# clump distance used
# get clump_dist_m used to create clumps in st_clump_points()
clump_dist_m <- st_clump_points_ans %>%
dplyr::mutate(clump_dist_m = as.numeric(clump_dist_m)) %>%
dplyr::filter(!is.na(clump_dist_m)) %>%
dplyr::slice(1) %>%
dplyr::pull(clump_dist_m)
if(length(clump_dist_m)!=1){stop("could not identify `clump_dist_m`")}
# get lb and ub levels
levels_clump_n_trees_lb <- levels(st_clump_points_ans$clump_n_trees_lb)
levels_clump_n_trees_ub <- levels(st_clump_points_ans$clump_n_trees_ub)
levels_clump_n_trees_grp <- levels(st_clump_points_ans$clump_n_trees_grp)
# numeric
num_clump_n_trees_lb <- levels_clump_n_trees_lb %>% as.character() %>% as.numeric()
num_clump_n_trees_ub <- levels_clump_n_trees_ub %>% as.character() %>% as.numeric()
##### we need to id the cut_breaks used in st_clump_points()
cut_breaks <-
c(
0 # always pad with 0
, num_clump_n_trees_ub %>%
unique() %>%
purrr::discard(~ is.null(.x) || is.na(.x))
, Inf # always pad with Inf
) %>%
unique() %>%
sort()
if(identical(cut_breaks,c(0,Inf))){
stop("could not identify cut breaks based on 'clump_n_trees_ub' factor levels")
}
return(list(
levels_clump_n_trees_lb = levels_clump_n_trees_lb
, levels_clump_n_trees_ub = levels_clump_n_trees_ub
, levels_clump_n_trees_grp = levels_clump_n_trees_grp
, num_clump_n_trees_lb = num_clump_n_trees_lb
, num_clump_n_trees_ub = num_clump_n_trees_ub
, cut_breaks = cut_breaks
, clump_dist_m = clump_dist_m
))
}
# check_clump_breaks(treetops_clustered)
cut_clump_fn <- function(
c # clump id from need_cut_trees data
, need_cut_trees # data with clumps already defined returned by st_clump_points()
, tgt # based on:
# as.numeric(as.character(levels(treetops_clustered$clump_n_trees_lb))) from st_clump_points()
# -or- min_clump_n_trees from get_target_check_prescription()
){
#### check target
tgt <- as.numeric(as.character(tgt))
if(!inherits(tgt,"numeric") || is.na(tgt)){
stop("`tgt` does not meet numeric conditions")
}
# does this exist as level in the data?
check_clump_breaks_ans <- check_clump_breaks(need_cut_trees)
lvl_exists <-
check_clump_breaks_ans$num_clump_n_trees_lb %>%
purrr::has_element(tgt)
if(!lvl_exists){
stop("`tgt` does not exist in 'clump_n_trees_lb' factor of `need_cut_trees` data")
}
##### we need to id the cut_breaks used in st_clump_points() to generate data in `need_cut_trees`
these_cut_breaks <- check_clump_breaks_ans$cut_breaks
# get distance used to create clumps in st_clump_points()
dist_temp <- check_clump_breaks_ans$clump_dist_m
# use only the current clump
curr_dta <- need_cut_trees %>%
dplyr::filter(clump_id == c)
if(nrow(curr_dta)<1){
stop("cannot find data with the clump_id == `c` in `need_cut_trees` data")
}
# do we even need to cut?
curr_clump_n_trees_lb <- curr_dta$clump_n_trees_lb %>%
as.character() %>%
as.numeric() %>%
unique()
if(length(curr_clump_n_trees_lb)!=1){
stop("multiple clumps identified in `need_cut_trees` data based on 'clump_n_trees_lb'")
}
# do we even need to cut?
if(
curr_clump_n_trees_lb<=tgt
|| nrow(curr_dta)<=tgt
){
message("no need to cut as clump contains fewer trees than `tgt`")
return(
curr_dta %>%
sf::st_drop_geometry() %>%
dplyr::mutate(is_keep_tree=1) %>%
dplyr::select(treeID, is_keep_tree)
)
}
# create clump polygon
clumps <- get_clump_summary(curr_dta)
# get the farthest line between points
f_line_temp <- st_points_to_line(curr_dta, line_ext = dist_temp)
# get perpendicular lines in dataset which we can iterate over to make cuts
perp_line_sf_temp <-
# for every 1 m along line length, get a new perp line
seq(
from = 0
, to = sf::st_length(f_line_temp) %>% as.numeric() %>% floor()
, by = 1
) %>%
purrr::map(
\(x)
st_perp_line(
interval=x, my_line = f_line_temp, proportion=F
)
) %>%
dplyr::bind_rows() %>%
dplyr::mutate(line_n = dplyr::row_number())
# find intersection of lines with the polygon and add length of intersection to perp line data
perp_line_sf_temp <- perp_line_sf_temp %>%
dplyr::inner_join(
# intersect and calc len
perp_line_sf_temp %>%
sf::st_intersection(
clumps %>%
dplyr::ungroup() %>%
dplyr::select(clump_id) %>%
dplyr::filter(clump_id == c)
) %>%
dplyr::mutate(len_m = sf::st_length(geometry) %>% as.numeric()) %>%
sf::st_drop_geometry()
, by = "line_n"
)
# make cuts at the points where there is the least overlap with the clump polygon
# list of potential line cut + tree combinations
# aggregate the mean length of the intersecting cut lines to the tree level
# prioritize trees for removal that have smallest length
cut_tree_lines_temp <-
curr_dta %>%
sf::st_buffer(dist_temp/2) %>%
sf::st_join(perp_line_sf_temp) %>%
sf::st_drop_geometry() %>%
dplyr::group_by(treeID) %>%
dplyr::summarise(len_m = mean(len_m, na.rm = T)) %>%
dplyr::ungroup() %>%
dplyr::arrange(len_m, treeID) %>%
dplyr::mutate(n = dplyr::row_number())
###############################################################
# while
###############################################################
# cut until the desired clump size is achieved based on these cut lines
while_temp <- 1
i <- 1
while(while_temp==1) {
# cut trees
cut_trees_temp <- cut_tree_lines_temp %>%
dplyr::slice(1:i) %>%
dplyr::distinct(treeID)
# get the remaining trees not cut
trees_remain_temp <- curr_dta %>%
dplyr::anti_join(cut_trees_temp, by = "treeID")
# ensure that there are trees
if(nrow(trees_remain_temp)==0){
if(best_desired_grps_n_temp==0){
# there are no more possible cuts :/
# ... so we're going to add trees until the desired clump size is reached
# start with biggest tree and keep adding trees until the desired clump size is reached
if( any(stringr::str_equal(names(curr_dta), "dbh_cm")) ){
start_tree <- curr_dta %>%
dplyr::arrange(desc(dbh_cm)) %>%
dplyr::slice(1)
}else if( any(stringr::str_equal(names(curr_dta), "tree_height_m")) ){
start_tree <- curr_dta %>%
dplyr::arrange(desc(tree_height_m)) %>%
dplyr::slice(1)
}else{
start_tree <- curr_dta %>% dplyr::slice(1)
}
increment_m <- 0.5
k <- 1
keep_trees <- dplyr::tibble(treeID = character(0))
while(nrow(keep_trees)==0){
# create a polygon to intersect with tree points and keep trees until number of trees met
poly_keep <- start_tree %>%
sf::st_buffer( (dist_temp/2) + increment_m*k ) %>%
dplyr::mutate(dummy = 1) %>%
dplyr::select(dummy)
# intersect and clump
i_trees <- curr_dta %>%
sf::st_intersection(poly_keep) %>%
st_clump_points(clump_dist_m = dist_temp, clump_breaks = these_cut_breaks)
# huh?
# i_trees %>% sf::st_drop_geometry() %>% dplyr::count(clump_n_trees_grp, clump_id)
# keep only the desired clump
keep_trees <- i_trees %>%
sf::st_drop_geometry() %>%
dplyr::filter(
clump_id ==
i_trees %>%
dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) == tgt) %>%
dplyr::pull(clump_id) %>%
.[1] %>%
dplyr::coalesce("nope")
) %>%
dplyr::select(treeID)
# increment
k <- k+1
}
# the remaining trees are cuts
best_cuts <- curr_dta %>%
sf::st_drop_geometry() %>%
dplyr::anti_join(keep_trees, by = "treeID") %>%
dplyr::select(treeID)
}else if(
best_desired_grps_n_temp>0
){
###############################
# add trees back in until gets worse
###############################
j <- 1
while_add <- 1
while(while_add==1){
# cut trees
cut_trees_temp = best_cuts %>%
# add trees back in (i.e. remove them from the cut trees)
dplyr::anti_join(
cut_tree_lines_temp %>%
dplyr::slice(1:j) %>%
dplyr::distinct(treeID)
, by = "treeID"
)
# get the remaining trees not cut
trees_remain_temp <- curr_dta %>%
dplyr::anti_join(cut_trees_temp, by = "treeID")
# count the groups remaining after cuts
desired_grps_n_temp <- trees_remain_temp %>%
st_clump_points(clump_dist_m = dist_temp, clump_breaks = these_cut_breaks) %>%
sf::st_drop_geometry() %>%
# do we have group sizes we want?
dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) == tgt) %>%
nrow()
if(desired_grps_n_temp>=best_desired_grps_n_temp){ # is this better than the best
best_cuts <- cut_trees_temp
best_desired_grps_n_temp <- desired_grps_n_temp
}else{
# stop it
while_add <- 0
}
# increment
j <- j + 1
} # while(while_add==1)
} # best_desired_grps_n_temp>0
# done so stop the whole stop it
while_temp <- 0
}else{ # if(nrow(trees_remain_temp)==0)
# count the groups remaining after cuts
desired_grps_n_temp <- trees_remain_temp %>%
st_clump_points(clump_dist_m = dist_temp, clump_breaks = these_cut_breaks) %>%
sf::st_drop_geometry() %>%
# do we have group sizes we want?
dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) == tgt) %>%
nrow()
### store best cut list
if(i==1){
best_cuts <- cut_trees_temp
best_desired_grps_n_temp <- desired_grps_n_temp
}else if(desired_grps_n_temp>best_desired_grps_n_temp){ # is this better than the best
best_cuts <- cut_trees_temp
best_desired_grps_n_temp <- desired_grps_n_temp
}else if(
desired_grps_n_temp==best_desired_grps_n_temp
&& i!=nrow(cut_tree_lines_temp)
){
best_cuts <- best_cuts
best_desired_grps_n_temp <- best_desired_grps_n_temp
}else if(
best_desired_grps_n_temp>0
&& desired_grps_n_temp<best_desired_grps_n_temp
){ # is this worse than the best which was successful?
###############################
# add trees back in until gets worse
###############################
j <- 1
while_add <- 1
while(while_add==1){
# cut trees
cut_trees_temp <- best_cuts %>%
# add trees back in (i.e. remove them from the cut trees)
dplyr::anti_join(
cut_tree_lines_temp %>%
dplyr::slice(1:j) %>%
dplyr::distinct(treeID)
, by = "treeID"
)
# get the remaining trees not cut
trees_remain_temp <- curr_dta %>%
dplyr::anti_join(cut_trees_temp, by = "treeID")
# count the groups remaining after cuts
desired_grps_n_temp <- trees_remain_temp %>%
st_clump_points(clump_dist_m = dist_temp, clump_breaks = these_cut_breaks) %>%
sf::st_drop_geometry() %>%
# do we have group sizes we want?
dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) == tgt) %>%
nrow()
if(desired_grps_n_temp>=best_desired_grps_n_temp){ # is this better than the best
best_cuts <- cut_trees_temp
best_desired_grps_n_temp <- desired_grps_n_temp
}else{
# stop it
while_add <- 0
}
# increment
j <- j + 1
} # while(while_add==1)
# done so stop the whole stop it
while_temp <- 0
}else if( i==nrow(cut_tree_lines_temp) ){ # is this the end?
# stop it
while_temp <- 0
}
### increment
i <- i+1
} # else if(nrow(trees_remain_temp)==0)
} # while(while_temp==1)
# return it
# return treelist with cut/keep
# join to original data and pull
d_temp <- curr_dta %>%
sf::st_drop_geometry() %>%
dplyr::mutate(is_keep_tree = as.numeric(NA)) %>%
dplyr::select(-c(is_keep_tree)) %>%
dplyr::left_join(
best_cuts %>% dplyr::mutate(is_keep_tree = 0)
, by = dplyr::join_by("treeID")
) %>%
dplyr::mutate(is_keep_tree = dplyr::coalesce(is_keep_tree, 1)) %>%
dplyr::select(treeID, is_keep_tree)
# returns treeID and keep tree flag data frame
return(d_temp)
}
# # example
# # pass it a clump id in the data
# check_clump_breaks_temp <- check_clump_breaks(treetops_clustered)
# # check_clump_breaks_temp
# xxx_temp <- treetops_clustered %>%
# dplyr::filter( as.numeric(as.character(clump_n_trees_lb)) == check_clump_breaks_temp$num_clump_n_trees_lb[3] ) %>%
# dplyr::arrange(desc(clump_n_trees)) %>%
# dplyr::slice(1) %>%
# # dplyr::select(tidyselect::starts_with("clump_"))
# dplyr::pull(clump_id) %>%
# cut_clump_fn(
# need_cut_trees = treetops_clustered
# , tgt = check_clump_breaks_temp$num_clump_n_trees_lb[2]
# )
# dplyr::glimpse(xxx_temp)
# # now we have to look at the clumps
# zzz_temp <-
# st_clump_points(
# x = treetops_sf_with_dbh %>%
# dplyr::inner_join(xxx_temp %>% dplyr::filter(is_keep_tree==1), by = "treeID")
# , clump_dist_m = check_clump_breaks_temp$clump_dist_m
# # , clump_breaks = these_cut_breaks
# , clump_breaks = check_clump_breaks_temp$cut_breaks
# )
# zzz_temp %>%
# dplyr::select(
# treeID, tidyselect::starts_with("clump_")
# ) %>%
# dplyr::glimpse()
# zzz_temp %>%
# # rename the active geometry column to "geometry"
# sf::st_set_geometry("geometry") %>%
# dplyr::mutate(geometry = sf::st_buffer(geometry, dist = clump_dist_m/2)) %>%
# ggplot2::ggplot() +
# ggplot2::geom_sf(
# mapping = ggplot2::aes(fill = clump_id)
# , alpha = 0.8
# ) +
# ggplot2::theme_void()there are some cases where the cut_clumps_fn may not
return the desired tree clumps based on a single pass using the initial
cut lines…iterate over the function until all clumps are in the desired
clump size or smaller.
get_cut_clump <- function(
c # clump id from need_cut_trees data
, need_cut_trees # data with clumps already defined returned by st_clump_points()
, tgt # based on:
# as.numeric(as.character(levels(treetops_clustered$clump_n_trees_lb))) from st_clump_points()
# -or- min_clump_n_trees from get_target_check_prescription()
) {
### check
check_clump_breaks_ans <- check_clump_breaks(need_cut_trees)
# we need to id the cut_breaks used in st_clump_points() to generate data in `need_cut_trees`
my_cut_breaks <- check_clump_breaks_ans$cut_breaks
# get distance used to create clumps in st_clump_points()
my_clump_dist_m <- check_clump_breaks_ans$clump_dist_m
#### check target
tgt <- as.numeric(as.character(tgt))
if(!inherits(tgt,"numeric") || is.na(tgt)){
stop("`tgt` does not meet numeric conditions")
}
# does this exist as level in the data?
lvl_exists <-
check_clump_breaks_ans$num_clump_n_trees_lb %>%
purrr::has_element(tgt)
if(!lvl_exists){
stop("`tgt` does not exist in 'clump_n_trees_lb' factor of `need_cut_trees` data")
}
# what if we need to keep cutting because could not find a solution based on initial cut lines?
cuts_first <- cut_clump_fn(
c = c
, need_cut_trees = need_cut_trees
, tgt = tgt
)
# did we get the target?
new_clumps <-
need_cut_trees %>%
dplyr::inner_join(
cuts_first %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
, by = "treeID"
) %>%
st_clump_points(clump_dist_m = my_clump_dist_m, clump_breaks = my_cut_breaks)
# new_clumps %>% sf::st_drop_geometry() %>% dplyr::count(clump_n_trees_grp)
while(
(
new_clumps %>% dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) > tgt) %>% nrow()
) > 0
){
# redo the cut clump
# redo the cut clump
cuts_again <- new_clumps %>%
dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) > tgt) %>%
dplyr::pull(clump_id) %>%
unique() %>%
purrr::map(
\(x)
cut_clump_fn(c = x, need_cut_trees = new_clumps, tgt = tgt)
) %>%
dplyr::bind_rows() %>%
dplyr::rename(updt = is_keep_tree)
# update the original cut data
cuts_first <- cuts_first %>%
dplyr::left_join(
cuts_again
, by = "treeID"
) %>%
dplyr::mutate(is_keep_tree = dplyr::coalesce(updt, is_keep_tree)) %>%
dplyr::select(-c(updt))
# reset the new clumps
new_clumps <-
need_cut_trees %>%
dplyr::inner_join(
cuts_first %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
, by = "treeID"
) %>%
st_clump_points(clump_dist_m = my_clump_dist_m, clump_breaks = my_cut_breaks)
# new_clumps %>% sf::st_drop_geometry() %>% dplyr::count(clump_n_trees_grp)
} # end while
return(cuts_first)
}function to pass a tree list
function to pass a whole tree list of sf point data
##############################################
# get cut keep tree flag
##############################################
get_keep_tree_flag <- function(
x # x = sf point data
, tgt # based on:
# as.numeric(as.character(levels(treetops_clustered$clump_n_trees_lb))) from st_clump_points()
# -or- min_clump_n_trees from get_target_check_prescription()
, clump_dist_m = NA # size (radius) of the epsilon neighborhood = maximum distance between points to add to clump
, clump_breaks = NA # where to break the clump size groups
, clump_labels = NA # what to name the clump size groups
) {
# MAKE CUTS INDEPENDENT OF CURRENT CLUMP GROUPING...
# what if we pass a tree list with trees already cut? and thus, potentially different clump group sizes than is
# currently defined in clump_n_trees_grp?
# 1) apply the clump grouping
# 2) only go through the cut alg for clumps that need cutting to the target
# 3) append the trees in target clump or lower to the list at the end
#### check target
tgt <- as.numeric(as.character(tgt))
if(!inherits(tgt,"numeric") || is.na(tgt)){
stop("`tgt` does not meet numeric conditions")
}
#############################
# 1) apply the clump grouping
#############################
new_clump_trees <- st_clump_points(
x = x
, clump_dist_m = clump_dist_m
, clump_breaks = clump_breaks
, clump_labels = clump_labels
# , ostory_ht_m = ostory_ht_m
# , ostory_dbh_cm = ostory_dbh_cm
)
### check
check_clump_breaks_ans <- check_clump_breaks(new_clump_trees)
# we need to id the cut_breaks used in st_clump_points() to generate data in `new_clump_trees`
my_cut_breaks <- check_clump_breaks_ans$cut_breaks
# get distance used to create clumps in st_clump_points()
my_clump_dist_m <- check_clump_breaks_ans$clump_dist_m
#### check target
# does this exist as level in the data?
lvl_exists <-
check_clump_breaks_ans$num_clump_n_trees_lb %>%
purrr::has_element(tgt)
if(!lvl_exists){
stop("`tgt` does not exist in 'clump_n_trees_lb' factor of `need_cut_trees` data")
}
# need cutting still
need_cut_trees <- new_clump_trees %>%
dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) > tgt)
#############################
# 2) only go through the cut alg for clumps that need cutting to the target
#############################
# for each clump that still needs cutting, iterate over and return keep/cut flag
compl_need_cut_trees <- need_cut_trees %>%
dplyr::pull(clump_id) %>%
unique() %>%
purrr::map(
\(x)
get_cut_clump(
c = x
, need_cut_trees = need_cut_trees
, tgt = tgt
)
) %>%
dplyr::bind_rows()
#############################
# 4) append the trees in target clump or lower to the list at the end
#############################
compl_need_cut_trees <- compl_need_cut_trees %>%
dplyr::bind_rows(
new_clump_trees %>%
sf::st_drop_geometry() %>%
dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) <= tgt) %>%
dplyr::mutate(is_keep_tree = 1) %>%
dplyr::select(treeID, is_keep_tree)
)
# return
return(
# original data so that order is preserved
x %>%
# add on trees that got cut with trees already good
dplyr::inner_join(
compl_need_cut_trees
, by = "treeID"
) %>%
dplyr::pull(is_keep_tree)
)
} # get_keep_tree_flagImplement Cutting
This is currently very messy and needs to be compiled into an easily accessible function…;/
cut groups incrementally to achieve target tree clump group size proportions
must first have:
- clustered trees from
get_tree_clumps()named astreetops_clustered - clump polygons from
get_clump_summary()named asclump_polys - target data from
get_target_check_prescription()named asprescription_target_data
# must first have:
#
# * clustered trees from `get_tree_clumps()`
# * clump polygons from `get_clump_summary()`
# * target data from `get_target_check_prescription()`
dir_temp <- "../data/"
fnm_temp <- file.path(dir_temp,"prescription_trees.gpkg") # what should the compiled rgb be called?
if(!dir.exists(dir_temp)){dir.create(dir_temp, showWarnings = F)}
if(!file.exists(fnm_temp)){
# check_clump_breaks(clump_polys) # will be the same
check_clump_breaks_temp <- check_clump_breaks(treetops_clustered) # will be the same
last_clump_n_trees_lb_temp <- min(check_clump_breaks_temp$num_clump_n_trees_lb)
# what is the smallest group with a target?
sm_grp_temp <-
prescription_target_data %>%
dplyr::filter(pct_stand_n_trees_target>0) %>%
dplyr::pull(min_clump_n_trees) %>%
# as.numeric(as.character(clump_n_trees_lb))
min()
# start with the largest group in target data currently with trees
grp_temp <-
prescription_target_data %>%
dplyr::arrange(desc(min_clump_n_trees)) %>%
dplyr::filter(pct_stand_n_trees_current>0) %>%
dplyr::slice(1) %>%
dplyr::pull(min_clump_n_trees)
# sm_grp_temp
# grp_temp
message(
paste("started cutting", grp_temp, "at", Sys.time())
)
# identify clumps to leave untouched
# clump_polys = get_clump_summary(ttops_temp)
keep_clumps_temp <-
clump_polys %>%
sf::st_drop_geometry() %>%
dplyr::mutate(
min_clump_n_trees = as.numeric(as.character(clump_n_trees_lb))
) %>%
#keep only trees in current group
dplyr::inner_join(
prescription_target_data %>%
dplyr::filter(min_clump_n_trees == grp_temp) %>%
dplyr::select(
mean_clump_n_trees, min_clump_n_trees, max_clump_n_trees
, stand_n_clumps_target
)
, by = "min_clump_n_trees"
) %>%
# keep only clumps that meet criteria
dplyr::filter(
n_trees >= min_clump_n_trees
& n_trees <= max_clump_n_trees
) %>%
# keep clumps closest the mean number in target
dplyr::mutate(
pct_to_target = abs( (n_trees-mean_clump_n_trees)/mean_clump_n_trees )
) %>%
dplyr::arrange(pct_to_target, desc(qmd_cm), n_trees, clump_id) %>%
dplyr::filter(dplyr::row_number()<=stand_n_clumps_target) %>%
dplyr::select(clump_id)
# dplyr::glimpse(keep_clumps_temp)
# !!!!!!!!!!FIX: WHAT IF WE HAVE CLUMPS OF THIS SIZE BUT THE N_TREES>MAX_TREES AND NEED TO GET CLUMPS OF THIS SIZE?
# ... UPDATE CUTTING ALG TO CUT CLUMP DOWN TO MIN-MAX TREE RANGE FOR CLUMPS WITH INF UPPER LIMIT
# start building tree list with keep/cut flag
# ttops_temp = get_tree_clumps(
# my_suid = my_suid
# , clump_dist_m = clump_dist_m
# , ostory_dbh_cm = ostory_dbh_cm
# )
# build tree list
keep_trees <-
treetops_clustered %>%
dplyr::ungroup() %>%
sf::st_drop_geometry() %>%
dplyr::inner_join(keep_clumps_temp, by = "clump_id") %>%
dplyr::select(treeID) %>%
dplyr::mutate(is_keep_tree = 1)
# dplyr::glimpse(keep_trees)
# did we keep all of the clumps?
if(
nrow(keep_clumps_temp) ==
( clump_polys %>% dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) == grp_temp) %>% nrow() )
){
remaining_trees <- dplyr::tibble(
treeID = character(0)
, is_keep_tree = numeric(0)
)
}else{
# determine keep/cut for the remaining trees in that group type
remaining_trees <-
treetops_clustered %>%
dplyr::ungroup() %>%
# keep only trees in current grp
dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) == grp_temp) %>%
# remove keep trees
dplyr::anti_join(keep_trees, by = "treeID")
# apply the cutting algorithm
# get the flag
remaining_trees$is_keep_tree <- get_keep_tree_flag(
x = remaining_trees
, tgt = prescription_target_data %>%
dplyr::ungroup() %>%
dplyr::arrange(min_clump_n_trees) %>%
dplyr::mutate(l = dplyr::lag(min_clump_n_trees)) %>%
dplyr::filter(min_clump_n_trees == grp_temp) %>%
dplyr::pull(l)
, clump_dist_m = check_clump_breaks_temp$clump_dist_m
, clump_breaks = check_clump_breaks_temp$cut_breaks
, clump_labels = check_clump_breaks_temp$levels_clump_n_trees_grp
)
# select relevant columns
remaining_trees <- remaining_trees %>%
sf::st_drop_geometry() %>%
dplyr::ungroup() %>%
dplyr::select(treeID, is_keep_tree)
}
message(
paste("done cutting", grp_temp, "at", Sys.time())
)
# remaining_trees %>% dplyr::count(is_keep_tree)
###############################################
# now process to go on to the next groups
###############################################
# get the next group
grp_temp <-
prescription_target_data %>%
dplyr::filter(min_clump_n_trees<grp_temp) %>%
dplyr::arrange(desc(min_clump_n_trees)) %>%
dplyr::slice(1) %>%
dplyr::pull(min_clump_n_trees)
while(grp_temp>=sm_grp_temp && grp_temp != last_clump_n_trees_lb_temp){
message(
paste("started cutting", grp_temp, "at", Sys.time())
)
# identify clumps to leave untouched
# clump_polys = get_clump_summary(ttops_temp)
keep_clumps_temp <-
clump_polys %>%
sf::st_drop_geometry() %>%
dplyr::mutate(
min_clump_n_trees = as.numeric(as.character(clump_n_trees_lb))
) %>%
#keep only trees in current group
dplyr::inner_join(
prescription_target_data %>%
dplyr::filter(min_clump_n_trees == grp_temp) %>%
dplyr::select(
mean_clump_n_trees, min_clump_n_trees, max_clump_n_trees
, stand_n_clumps_target
)
, by = "min_clump_n_trees"
) %>%
# keep only clumps that meet criteria
dplyr::filter(
n_trees >= min_clump_n_trees
& n_trees <= max_clump_n_trees
) %>%
# keep clumps closest the mean number in target group
dplyr::mutate(
pct_to_target = abs( (n_trees-mean_clump_n_trees)/mean_clump_n_trees )
) %>%
dplyr::arrange(pct_to_target, desc(qmd_cm), n_trees, clump_id) %>%
dplyr::filter(dplyr::row_number()<=stand_n_clumps_target) %>%
dplyr::select(clump_id)
# add to tree list with keep/cut flag
keep_trees <- keep_trees %>%
dplyr::bind_rows(
treetops_clustered %>%
dplyr::ungroup() %>%
sf::st_drop_geometry() %>%
dplyr::inner_join(keep_clumps_temp, by = "clump_id") %>%
dplyr::select(treeID) %>%
dplyr::mutate(is_keep_tree = 1)
)
# check if the desired clump number was met and add the remaining trees from previous group if needed
more_clumps_target <- 0
if(
nrow(keep_clumps_temp) <
(
prescription_target_data %>%
dplyr::filter(min_clump_n_trees == grp_temp) %>%
dplyr::pull(stand_n_clumps_target)
)
){
# how many more clumps are needed?
more_clumps_target <-
(
prescription_target_data %>%
dplyr::filter(min_clump_n_trees == grp_temp) %>%
dplyr::pull(stand_n_clumps_target)
) - nrow(keep_clumps_temp)
# determine group size of remaining trees in the group prior that were not in a group selected and were not cut
potential_trees <-
treetops_clustered %>%
dplyr::ungroup() %>%
dplyr::inner_join(
remaining_trees %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
, by = "treeID"
) %>%
st_clump_points(
clump_dist_m = check_clump_breaks_temp$clump_dist_m
, clump_breaks = check_clump_breaks_temp$cut_breaks
, clump_labels = check_clump_breaks_temp$levels_clump_n_trees_grp
) %>%
# ggplot() + geom_sf(aes(fill = clump_n_trees_grp)) + theme_void()
# get the original clump id to prioritize new clump groups in the same area
dplyr::inner_join(
treetops_clustered %>%
sf::st_drop_geometry() %>%
dplyr::ungroup() %>%
dplyr::select(treeID, clump_id) %>%
dplyr::rename(orig_clump_id = clump_id)
, by = "treeID"
) %>%
dplyr::group_by(orig_clump_id) %>%
dplyr::mutate(
pct_desired_grp =
sum(ifelse( as.numeric(as.character(clump_n_trees_lb)) == grp_temp, 1, 0)) / dplyr::n()
) %>%
dplyr::ungroup() %>%
# keep only the current group
dplyr::filter( as.numeric(as.character(clump_n_trees_lb)) == grp_temp)
# pick trees from potential trees based on clump summary
keep_trees <- keep_trees %>%
dplyr::bind_rows(
potential_trees %>%
sf::st_drop_geometry() %>%
# filter trees based on clumps needed
dplyr::inner_join(
get_clump_summary(potential_trees) %>%
sf::st_drop_geometry() %>%
# join with original clump id metrics to prioritize
# keeping clumps in the same area and minimize cutting time
dplyr::left_join(
potential_trees %>%
sf::st_drop_geometry() %>%
dplyr::group_by(clump_id, orig_clump_id) %>%
dplyr::summarise(pct_desired_grp = max(pct_desired_grp))
, by = "clump_id"
) %>%
# keep the number of clumps needed
dplyr::arrange(desc(pct_desired_grp), orig_clump_id, desc(qmd_cm), n_trees, clump_id) %>%
dplyr::filter(dplyr::row_number()<=more_clumps_target) %>%
dplyr::select(clump_id)
, by = "clump_id"
) %>%
dplyr::select(treeID) %>%
dplyr::mutate(is_keep_tree = 1)
)
} # end if don't have enough clumps
####################################################
# determine keep/cut for the remaining trees in that group type and higher groups
####################################################
remaining_trees <-
treetops_clustered %>%
dplyr::ungroup() %>%
# REMOVE TREES FROM PREVIOUS TREES REMAINING that got cut to make this group size
dplyr::anti_join(
remaining_trees %>% dplyr::filter(is_keep_tree == 0)
, by = "treeID"
) %>%
# keep only trees in current grp or prior grp
dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) >= grp_temp) %>%
# remove keep trees
dplyr::anti_join(keep_trees, by = "treeID")
if(nrow(remaining_trees) > 0){
# apply the cutting algorithm
# get the flag
remaining_trees$is_keep_tree <- get_keep_tree_flag(
x = remaining_trees
, tgt = prescription_target_data %>%
dplyr::ungroup() %>%
dplyr::arrange(min_clump_n_trees) %>%
dplyr::mutate(l = dplyr::lag(min_clump_n_trees)) %>%
dplyr::filter(min_clump_n_trees == grp_temp) %>%
dplyr::pull(l)
, clump_dist_m = check_clump_breaks_temp$clump_dist_m
, clump_breaks = check_clump_breaks_temp$cut_breaks
, clump_labels = check_clump_breaks_temp$levels_clump_n_trees_grp
)
# select relevant columns
remaining_trees <- remaining_trees %>%
sf::st_drop_geometry() %>%
dplyr::ungroup() %>%
dplyr::select(treeID, is_keep_tree)
}
# increment
message(
paste("done cutting", grp_temp, "at", Sys.time())
)
# get the next group
grp_temp <-
prescription_target_data %>%
dplyr::filter(min_clump_n_trees<grp_temp) %>%
dplyr::arrange(desc(min_clump_n_trees)) %>%
dplyr::slice(1) %>%
dplyr::pull(min_clump_n_trees) %>%
dplyr::coalesce(last_clump_n_trees_lb_temp)
}
# now individuals
if(sm_grp_temp == last_clump_n_trees_lb_temp){
message(
paste("started cutting", sm_grp_temp, "at", Sys.time())
)
potential_trees <-
# original data
treetops_clustered %>%
dplyr::filter( as.numeric(as.character(clump_n_trees_lb)) == last_clump_n_trees_lb_temp) %>%
# remaining trees
dplyr::bind_rows(
treetops_clustered %>%
dplyr::ungroup() %>%
dplyr::inner_join(
remaining_trees %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
, by = "treeID"
)
) %>%
# make sure that these are all individuals
dplyr::mutate(is_orig = ifelse( as.numeric(as.character(clump_n_trees_lb))==last_clump_n_trees_lb_temp, 1, 0)) %>%
st_clump_points(
clump_dist_m = check_clump_breaks_temp$clump_dist_m
, clump_breaks = check_clump_breaks_temp$cut_breaks
, clump_labels = check_clump_breaks_temp$levels_clump_n_trees_grp
) %>%
dplyr::filter( as.numeric(as.character(clump_n_trees_lb)) == last_clump_n_trees_lb_temp)
# pick trees from potential trees based on target
keep_trees <- keep_trees %>%
dplyr::select(treeID) %>%
dplyr::bind_rows(
potential_trees %>%
sf::st_drop_geometry() %>%
# keep the number of clumps needed
dplyr::arrange(desc(is_orig), desc(dbh_cm), desc(tree_height_m)) %>%
dplyr::filter(
dplyr::row_number() <=
(
prescription_target_data %>%
dplyr::filter(min_clump_n_trees == last_clump_n_trees_lb_temp) %>%
dplyr::pull(stand_n_clumps_target)
)
) %>%
dplyr::select(treeID)
) %>%
dplyr::mutate(is_keep_tree = 1)
message(
paste("done cutting", sm_grp_temp, "at", Sys.time())
)
}
# get the final prescription
# treetops_clustered %>% dplyr::select(tidyselect::starts_with("clump_n_trees_")) %>% names()
# treetops_clustered %>%
# dplyr::rename_with(
# .cols = tidyselect::starts_with("clump_")
# , .fn = ~ paste0("orig_", .x, recycle0 = TRUE)
# ) %>%
# names()
prescription_trees <-
treetops_clustered %>%
dplyr::ungroup() %>%
dplyr::left_join(
keep_trees
, by = "treeID"
) %>%
# tracking vars
dplyr::mutate(
# fill in keep tree flag
is_keep_tree = dplyr::coalesce(is_keep_tree, 0)
) %>%
dplyr::rename_with(
.cols = tidyselect::starts_with("clump_")
, .fn = ~ paste0("orig_", .x, recycle0 = TRUE)
)
# attach the new clumping to the trees
prescription_trees <- prescription_trees %>%
dplyr::left_join(
# reclump
prescription_trees %>%
dplyr::filter(is_keep_tree == 1) %>%
st_clump_points(
clump_dist_m = check_clump_breaks_temp$clump_dist_m
, clump_breaks = check_clump_breaks_temp$cut_breaks
, clump_labels = check_clump_breaks_temp$levels_clump_n_trees_grp
) %>%
sf::st_drop_geometry() %>%
dplyr::select(treeID, tidyselect::starts_with("clump_")) %>%
dplyr::rename_with(
.cols = tidyselect::starts_with("clump_")
, .fn = ~ paste0("new_", .x, recycle0 = TRUE)
)
, by = "treeID"
) %>%
dplyr::mutate(
new_clump_n_trees_grp = forcats::fct_na_value_to_level(new_clump_n_trees_grp, level = "Cut tree")
)
# save
prescription_trees %>% sf::st_write(fnm_temp, append = F)
}else{
prescription_trees <- sf::st_read(fnm_temp, quiet = T)
}check the distribution of current vs prescription tree clump sizes
prescription_trees %>%
sf::st_drop_geometry() %>%
dplyr::count(
orig_clump_n_trees_grp, new_clump_n_trees_grp
) %>%
tidyr::pivot_wider(names_from = orig_clump_n_trees_grp, values_from = n, values_fill = 0) %>%
dplyr::arrange(new_clump_n_trees_grp) %>%
dplyr::rename(`New` = new_clump_n_trees_grp) %>%
kableExtra::kbl() %>%
kableExtra::kable_styling() %>%
kableExtra::add_header_above(
c(" "=1,"Old"=length(unique(prescription_trees$orig_clump_n_trees_grp)))
)| New | 10-15 trees | 16-25 trees | 2-4 trees | 5-9 trees | >25 trees | Individual |
|---|---|---|---|---|---|---|
| 10-15 trees | 356 | 65 | 0 | 0 | 325 | 0 |
| 16-25 trees | 0 | 35 | 0 | 0 | 139 | 0 |
| 2-4 trees | 0 | 70 | 731 | 0 | 701 | 0 |
| 5-9 trees | 0 | 38 | 0 | 607 | 514 | 0 |
| >25 trees | 0 | 0 | 0 | 0 | 228 | 0 |
| Cut tree | 0 | 184 | 0 | 0 | 3844 | 0 |
| Individual | 0 | 22 | 0 | 0 | 595 | 354 |
get the distance raster and openings vector data
treetops_clustered_new <-
prescription_trees %>%
dplyr::filter(is_keep_tree == 1) %>%
dplyr::select(-tidyselect::starts_with("orig_")) %>%
# sf::st_drop_geometry() %>% dplyr::count(new_clump_n_trees_grp)
dplyr::rename_with(
.cols = tidyselect::starts_with("new_clump_")
, .fn = ~ stringr::str_remove(.x, "^new_")
)
# dplyr::glimpse(treetops_clustered_new)
clump_polys_new <- get_clump_summary(treetops_clustered_new)
# dplyr::glimpse(clump_polys_new)
# clump_polys_new %>% dplyr::select(clump_id) %>% plot()
get_clump_dist_rast_ans2 <- get_clump_dist_rast(get_clump_summary_ans = clump_polys_new, stand_sf = stand_sf)check out the openings information
# what are the openings
get_clump_dist_rast_ans2$openings_vect$openining_area_m2 %>% summary()## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.04 0.04 1.92 71.53 17.96 11077.28get_clump_dist_rast_ans2$dist_rast %>% terra::aggregate(2, cores = 4) %>% terra::summary()## |---------|---------|---------|---------|========================================= ## layer
## Min. : 0.000
## 1st Qu.: 0.000
## Median : 0.714
## Mean : 1.641
## 3rd Qu.: 2.500
## Max. :21.877
## NA's :29973plot the prescription
plt_fnl_12 <-
# ggplot2::ggplot() +
plt_stand +
# distance
ggplot2::geom_tile(
data = get_clump_dist_rast_ans2$dist_rast %>% terra::aggregate(2, cores = 4) %>% terra::as.data.frame(xy = T) %>% dplyr::rename(f=3)
, mapping = ggplot2::aes(x=x, y=y, fill = f)
) +
harrypotter::scale_fill_hp(
"mischief"
, na.value = "transparent"
, name = "distance to\nnearest tree (m)"
, limits = c(0,22)
) +
# openings
ggplot2::geom_sf(
data = get_clump_dist_rast_ans2$openings_vect
, mapping = ggplot2::aes(color = openining_area_m2), fill = NA
) +
ggplot2::scale_color_gradient(
low = "gray77", high = "gray11"
, labels = scales::comma_format(accuracy = 1)
, name = latex2exp::TeX("opening\narea ($\\m^2$)")
, limits = c(0,11100)
) +
# clumps
ggnewscale::new_scale_fill() +
ggplot2::geom_sf(
data = clump_polys_new
, mapping = ggplot2::aes(fill = clump_n_trees_grp)
, color = NA
) +
harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
# tree points
ggplot2::geom_sf(data = treetops_clustered_new, color = "gray88", shape = ".") +
ggplot2::theme_void() +
ggplot2::theme(
legend.position = "right"
, legend.direction = "vertical"
, plot.subtitle = ggplot2::element_text(size = 10, hjust = 0.5, face = "bold")
, legend.title= ggplot2::element_text(size=8)
, legend.text= ggplot2::element_text(size = 7)
)## |---------|---------|---------|---------|========================================= plot it
# plot
plt_fnl_12
# save it
ggplot2::ggsave("../data/06_clumps_dist_presciption.jpg", width = 7, height = 5)highlight the openings
plt_open_12 <-
plt_stand +
# clumps
ggplot2::geom_sf(
data = clump_polys_new
, mapping = ggplot2::aes(fill = clump_n_trees_grp)
, color = NA
) +
harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
# openings
ggnewscale::new_scale_fill() +
ggplot2::geom_sf(
data = get_clump_dist_rast_ans2$openings_vect
, mapping = ggplot2::aes(fill = openining_area_m2), color = NA
) +
ggplot2::scale_fill_gradient(
low = "gray77", high = "gray11"
, labels = scales::comma_format(accuracy = 1)
, name = latex2exp::TeX("opening\narea ($\\m^2$)")
, limits = c(0,11100)
) +
# tree points
ggplot2::geom_sf(data = treetops_clustered_new, color = "gray88", shape = ".") +
ggplot2::theme_void() +
ggplot2::theme(
legend.position = "right"
, legend.direction = "vertical"
, plot.subtitle = ggplot2::element_text(size = 10, hjust = 0.5, face = "bold")
, legend.title= ggplot2::element_text(size=8)
, legend.text= ggplot2::element_text(size = 7)
)
# plot
plt_open_12
ggplot2::ggsave("../data/07_clumps_openings_presciption.jpg", width = 7, height = 5)combine them with patchwork
plt_fnl_12 + (plt_open_12 + theme(legend.position = "none")) 
# &
# ggplot2::theme(
# plot.subtitle = ggplot2::element_text(size = 10, hjust = 0.5, face = "bold")
# , legend.title=ggplot2::element_text(size=8)
# , legend.text=ggplot2::element_text(size = 7)
# )
ggplot2::ggsave("../data/08_clumps_combined_prescription.jpg", width = 10, height = 5.5)combine old vs new in plot
(plt_fnl_1 + labs(subtitle = "Pre-Treatment")) +
(plt_fnl_12 + labs(subtitle = "Post-Treatment") + theme(legend.position = "none")) &
theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))
# save it
ggplot2::ggsave("../data/09_clumps_dist_compared.jpg", width = 10, height = 7)combine old vs new openings highlight
(plt_open_1 + labs(subtitle = "Pre-Treatment")) +
(plt_open_12 + labs(subtitle = "Post-Treatment") + theme(legend.position = "none")) &
theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))
ggplot2::ggsave("../data/10_clumps_openings_compared.jpg", width = 10, height = 7)table comparison
current vs target vs prescription
clump_n_trees_grp_summary2 <-
get_clump_n_trees_grp_summary(
get_tree_clumps_ans = treetops_clustered_new
, get_clump_summary_ans = clump_polys_new
, stand_sf = stand_sf
)
# huh?
clump_n_trees_grp_summary2 %>% dplyr::glimpse()## Rows: 6
## Columns: 33
## $ stand_area_ha <dbl> 30.4578, 30.4578, 30.4578, 30.4578, 30.457…
## $ clump_n_trees_grp <chr> "10-15 trees", "16-25 trees", "2-4 trees",…
## $ clump_n_trees_lb <chr> "10", "16", "2", "5", "26", "1"
## $ clump_n_trees_ub <chr> "15", "25", "4", "9", "Inf", "1"
## $ n_trees <int> 746, 174, 1502, 1159, 228, 971
## $ mean_dbh_cm <dbl> 27.14707, 28.91884, 27.77161, 27.56157, 27…
## $ mean_tree_height_m <dbl> 14.71029, 16.23725, 15.14332, 14.96687, 15…
## $ loreys_height_m <dbl> 15.96605, 17.23514, 16.21269, 16.05672, 16…
## $ basal_area_m2 <dbl> 44.93656, 11.72872, 93.86354, 71.45057, 14…
## $ clump_area_ha <dbl> 1.631763, 0.376619, 3.562791, 2.597137, 0.…
## $ stand_n_clumps <int> 63, 9, 513, 185, 7, 971
## $ trees_per_ha <dbl> 24.492902, 5.712822, 49.314127, 38.052645,…
## $ basal_area_m2_per_ha <dbl> 1.4753710, 0.3850810, 3.0817568, 2.3458871…
## $ qmd_cm <dbl> 27.69399, 29.29583, 28.20776, 28.01665, 28…
## $ stand_trees_per_ha <dbl> 156.9384, 156.9384, 156.9384, 156.9384, 15…
## $ stand_basal_area_m2 <dbl> 303.453, 303.453, 303.453, 303.453, 303.45…
## $ stand_basal_area_m2_per_ha <dbl> 9.963063, 9.963063, 9.963063, 9.963063, 9.…
## $ pct_stand_basal_area <dbl> 0.14808408, 0.03865086, 0.30931822, 0.2354…
## $ pct_stand_n_trees <dbl> 0.15606695, 0.03640167, 0.31422594, 0.2424…
## $ stand_qmd_cm <dbl> 28.43065, 28.43065, 28.43065, 28.43065, 28…
## $ mean_dbh_in <dbl> 10.69594, 11.39402, 10.94202, 10.85926, 11…
## $ qmd_in <dbl> 10.91143, 11.54256, 11.11386, 11.03856, 11…
## $ stand_qmd_in <dbl> 11.20168, 11.20168, 11.20168, 11.20168, 11…
## $ mean_tree_height_ft <dbl> 48.24976, 53.25817, 49.67010, 49.09134, 50…
## $ loreys_height_ft <dbl> 52.36863, 56.53127, 53.17761, 52.66603, 53…
## $ basal_area_ft2_per_ac <dbl> 6.431142, 1.678568, 13.433378, 10.225722, …
## $ stand_basal_area_ft2_per_ac <dbl> 43.42899, 43.42899, 43.42899, 43.42899, 43…
## $ trees_per_ac <dbl> 9.919625, 2.313693, 19.972221, 15.411321, …
## $ stand_trees_per_ac <dbl> 63.56006, 63.56006, 63.56006, 63.56006, 63…
## $ stand_area_ac <dbl> 75.26123, 75.26123, 75.26123, 75.26123, 75…
## $ clump_area_ac <dbl> 4.0320853, 0.9306256, 8.8036559, 6.4175263…
## $ basal_area_ft2 <dbl> 483.6971, 126.2479, 1010.3472, 769.0939, 1…
## $ stand_basal_area_ft2 <dbl> 3266.368, 3266.368, 3266.368, 3266.368, 32…table it to compare targets versus prescription
prescription_target_data %>%
dplyr::select(
clump_n_trees_grp
, tidyselect::starts_with("pct_stand_n_trees")
, tidyselect::starts_with("trees_per_acre_")
) %>%
dplyr::left_join(
clump_n_trees_grp_summary2 %>%
dplyr::select(clump_n_trees_grp, pct_stand_n_trees, trees_per_ac) %>%
dplyr::rename(
pct_stand_n_trees_presc = pct_stand_n_trees
, trees_per_acre_presc = trees_per_ac
)
, by = "clump_n_trees_grp"
) %>%
dplyr::mutate(
dplyr::across(
tidyselect::where(is.numeric)
, ~ dplyr::coalesce(.x,0)
)
) %>%
dplyr::relocate(
pct_stand_n_trees_presc
, .after = pct_stand_n_trees_target
) %>%
dplyr::relocate(
trees_per_acre_presc
, .after = trees_per_acre_target
) %>%
dplyr::mutate(
dplyr::across(
tidyselect::starts_with("pct_stand_n_trees")
, ~ scales::percent(.x,accuracy = 1)
)
) %>%
kableExtra::kable(
caption = paste0(
"Current stand (ft<sup>2</sup> ac<sup>-1</sup>): "
, clump_n_trees_grp_summary$stand_basal_area_ft2_per_ac[1] %>% scales::number(accuracy = 0.1)
, "<br>Current stand QMD (in): "
, clump_n_trees_grp_summary$stand_qmd_in[1] %>% scales::number(accuracy = 0.1)
, "<br>Current stand TPA: "
, clump_n_trees_grp_summary$stand_trees_per_ac[1] %>% scales::number(accuracy = 1)
)
, escape = F
, digits = 1
, col.names = c(
"", rep(c("current","target","prescription"),2)
)
) %>%
kableExtra::kable_styling() %>%
kableExtra::add_header_above(
c(" "=1,"% Trees"=3, "TPA"=3)
)| current | target | prescription | current | target | prescription | |
|---|---|---|---|---|---|---|
| Individual | 4% | 25% | 20% | 4.7 | 16.0 | 12.9 |
| 2-4 trees | 8% | 40% | 31% | 9.7 | 25.6 | 20.0 |
| 5-9 trees | 7% | 25% | 24% | 8.1 | 16.0 | 15.4 |
| 10-15 trees | 4% | 10% | 16% | 4.7 | 6.4 | 9.9 |
| 16-25 trees | 5% | 0% | 4% | 5.5 | 0.0 | 2.3 |
| >25 trees | 72% | 0% | 5% | 84.4 | 0.0 | 3.0 |
check the pre- and post-treatment summary
# table it
clump_n_trees_grp_summary2 %>%
dplyr::mutate(presc = "post-treatment") %>%
dplyr::bind_rows(
clump_n_trees_grp_summary %>%
dplyr::mutate(presc = "pre-treatment")
) %>%
dplyr::arrange(clump_n_trees_grp, desc(presc)) %>%
dplyr::select(
clump_n_trees_grp, presc
, n_trees
, mean_dbh_in
, qmd_in
, mean_tree_height_ft
, loreys_height_ft
, trees_per_ac
, basal_area_ft2_per_ac, pct_stand_basal_area, pct_stand_n_trees
) %>%
dplyr::mutate(
dplyr::across(
.cols = c(pct_stand_basal_area, pct_stand_n_trees)
, .fns = ~ scales::percent(.x, accuracy = 1)
)
) %>%
kableExtra::kbl(
digits = 1
, escape = F
, caption = "Pre- and Post-treatment overstory tree clump summary"
, col.names = c(
".", ""
, "trees"
, "mean<br>DBH (in)"
, "QMD (in)"
, "mean<br>Ht. (ft)"
, "Loreys<br>Ht. (ft)"
, "TPA"
, "BA<br>ft<sup>2</sup> ac<sup>-1</sup>"
, "%<br>stand BA"
, "%<br>stand trees"
)
) %>%
kableExtra::kable_styling() %>%
kableExtra::collapse_rows(columns = 1, valign = "top")| . | trees |
mean DBH (in) |
QMD (in) |
mean Ht. (ft) |
Loreys Ht. (ft) |
TPA |
BA ft2 ac-1 |
% stand BA |
% stand trees |
|
|---|---|---|---|---|---|---|---|---|---|---|
| 10-15 trees | pre-treatment | 356 | 10.5 | 10.8 | 47.7 | 52.8 | 4.7 | 3.0 | 4% | 4% |
| post-treatment | 746 | 10.7 | 10.9 | 48.2 | 52.4 | 9.9 | 6.4 | 15% | 16% | |
| 16-25 trees | pre-treatment | 414 | 10.5 | 10.7 | 46.9 | 50.8 | 5.5 | 3.5 | 4% | 5% |
| post-treatment | 174 | 11.4 | 11.5 | 53.3 | 56.5 | 2.3 | 1.7 | 4% | 4% | |
| 2-4 trees | pre-treatment | 731 | 10.9 | 11.1 | 49.8 | 54.1 | 9.7 | 6.5 | 8% | 8% |
| post-treatment | 1502 | 10.9 | 11.1 | 49.7 | 53.2 | 20.0 | 13.4 | 31% | 31% | |
| 5-9 trees | pre-treatment | 607 | 10.6 | 10.8 | 47.9 | 52.3 | 8.1 | 5.2 | 7% | 7% |
| post-treatment | 1159 | 10.9 | 11.0 | 49.1 | 52.7 | 15.4 | 10.2 | 24% | 24% | |
| >25 trees | pre-treatment | 6346 | 10.9 | 11.0 | 48.6 | 51.4 | 84.4 | 55.7 | 72% | 72% |
| post-treatment | 228 | 11.0 | 11.2 | 50.1 | 53.6 | 3.0 | 2.1 | 5% | 5% | |
| Individual | pre-treatment | 354 | 10.9 | 11.1 | 50.4 | 55.1 | 4.7 | 3.2 | 4% | 4% |
| post-treatment | 971 | 11.6 | 11.7 | 53.7 | 55.6 | 12.9 | 9.6 | 22% | 20% |